Le problème d’apprentissage

Le chapitre précédent a présenté les méthodes non paramétriques, comme les k plus proches voisins, qui conservent les données d’entraînement et les consultent au moment de la prédiction. Cette approche est intuitive, mais elle a un coût: les données doivent rester en mémoire, et chaque prédiction requiert de parcourir l’ensemble d’entraînement. Ce chapitre développe une approche différente: plutôt que de garder les données, nous cherchons à les résumer dans un ensemble fixe de paramètres. L’apprentissage devient alors un problème d’optimisation.

Apprentissage supervisé¶

Une ingénieure automobile mesure la distance de freinage d’un véhicule à différentes vitesses. Ses données ressemblent à ceci:

Elle veut prédire la distance de freinage à 30 mph sans faire l’essai. Pour cela, elle cherche une fonction $f$ telle que $f(\text{vitesse}) \approx \text{distance}$ sur ses observations. Si la fonction capture la relation sous-jacente, elle devrait donner une prédiction raisonnable pour des vitesses non mesurées.

Source

import numpy as np
import matplotlib.pyplot as plt

# Configuration pour des figures haute résolution
%config InlineBackend.figure_format = 'retina'

# Données de freinage (Ezekiel, 1930): vitesse (mph) vs distance d'arrêt (ft)
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

plt.figure(figsize=(6, 4))
plt.scatter(speed, dist, alpha=0.7, label='Observations')

# Fit quadratic
coeffs = np.polyfit(speed, dist, 2)
speed_grid = np.linspace(0, 30, 100)
dist_pred = np.polyval(coeffs, speed_grid)
plt.plot(speed_grid, dist_pred, 'k--', alpha=0.6, label='Fonction ajustée')

# Prediction at 30 mph
pred_30 = np.polyval(coeffs, 30)
plt.scatter([30], [pred_30], marker='x', s=80, color='C1', zorder=5, label=f'Prédiction à 30 mph: {pred_30:.0f} ft')

plt.xlabel('Vitesse (mph)')
plt.ylabel('Distance de freinage (ft)')
plt.legend()
plt.tight_layout()

Ce processus est l’ajustement de courbe (curve fitting). Nous avons des paires (entrée, sortie), nous ajustons une fonction, et nous utilisons cette fonction pour prédire. L’apprentissage supervisé généralise cette idée: les entrées peuvent être des vecteurs de dimension quelconque, les sorties peuvent être continues ou discrètes, et les fonctions candidates peuvent être bien plus complexes qu’un polynôme.

Formellement, nous disposons d’un jeu de données $\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^N$ composé de $N$ paires, où chaque $\mathbf{x}_i \in \mathcal{X}$ est une entrée et $y_i \in \mathcal{Y}$ est la sortie correspondante. L’objectif est de trouver une fonction $f: \mathcal{X} \to \mathcal{Y}$ qui approxime bien la relation entre entrées et sorties, y compris pour des exemples que nous n’avons pas encore observés.

Dans de nombreuses applications, les entrées sont des vecteurs de caractéristiques. Chaque exemple $\mathbf{x}_i \in \mathbb{R}^d$ est un vecteur de dimension $d$, où chaque composante représente une mesure ou un attribut. Pour prédire le prix d’une maison, les entrées pourraient être la superficie, le nombre de chambres et l’âge du bâtiment. Pour filtrer les pourriels, les entrées pourraient être des fréquences de mots. Pour diagnostiquer une maladie, les entrées pourraient être des résultats d’analyses sanguines.

Lorsque la sortie est une valeur continue, nous parlons de régression: $f: \mathbb{R}^d \to \mathbb{R}$ pour une sortie scalaire, ou $f: \mathbb{R}^d \to \mathbb{R}^p$ pour une sortie vectorielle. La distance de freinage, le prix d’une maison, la concentration d’un médicament dans le sang sont des exemples de régression.

Lorsque la sortie est une catégorie parmi un ensemble fini, nous parlons de classification. Pour la classification binaire, $f: \mathbb{R}^d \to \{0, 1\}$. Pour la classification multiclasse avec $m$ catégories, $f: \mathbb{R}^d \to \{0, \ldots, m-1\}$. Déterminer si un courriel est un pourriel, diagnostiquer une maladie, ou reconnaître un chiffre manuscrit sont des exemples de classification.

Mesurer l’erreur¶

Pour choisir entre deux fonctions candidates, nous avons besoin d’un critère qui quantifie la qualité des prédictions. Une fonction de perte $\ell: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}_+$ mesure l’écart entre une prédiction $\hat{y}$ et la vraie valeur $y$. Une perte de zéro indique une prédiction parfaite; plus la perte est grande, plus l’erreur est importante.

Pour la régression, nous utilisons généralement la perte quadratique:

Cette perte pénalise les grandes erreurs de manière quadratique. Une erreur de 2 coûte quatre fois plus qu’une erreur de 1.

Reprenons les données de freinage. Supposons que notre fonction prédise 50 ft pour une vitesse où la vraie distance est 56 ft. La perte quadratique est $(56 - 50)^2 = 36$. Si elle prédit 70 ft, la perte est $(56 - 70)^2 = 196$. La perte quadratique pénalise sévèrement les grandes erreurs.

Source

import numpy as np
import matplotlib.pyplot as plt

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

coeffs = np.polyfit(speed, dist, 2)
predictions = np.polyval(coeffs, speed)
residuals = dist - predictions

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# Left: predictions vs observations
ax = axes[0]
ax.scatter(speed, dist, alpha=0.7, label='Observations')
speed_grid = np.linspace(4, 25, 100)
ax.plot(speed_grid, np.polyval(coeffs, speed_grid), 'k--', alpha=0.6, label='Prédictions')
for i in range(len(speed)):
    ax.plot([speed[i], speed[i]], [dist[i], predictions[i]], 'C1-', alpha=0.5)
ax.set_xlabel('Vitesse (mph)')
ax.set_ylabel('Distance (ft)')
ax.legend()
ax.set_title('Résidus: écarts entre observations et prédictions')

# Right: histogram of squared residuals
ax = axes[1]
ax.hist(residuals**2, bins=15, edgecolor='black', alpha=0.7)
ax.set_xlabel(r'Perte quadratique $(y - \hat{y})^2$')
ax.set_ylabel('Fréquence')
ax.set_title(f'EQM = {np.mean(residuals**2):.1f}')

plt.tight_layout()

Pour la classification, un choix naturel est la perte 0-1:

Cette perte compte simplement les erreurs: elle vaut 1 pour une mauvaise prédiction, 0 sinon.

Le choix de la fonction de perte dépend du problème. En diagnostic médical, manquer une maladie grave (faux négatif) peut avoir des conséquences bien plus importantes que de prescrire un test supplémentaire à un patient sain (faux positif). Une perte asymétrique refléterait cette différence. En régression, si les grandes erreurs sont particulièrement problématiques, la perte quadratique est appropriée; si nous voulons être robustes aux valeurs aberrantes, la perte absolue $|y - \hat{y}|$ est préférable.

Le risque¶

La perte évalue une seule prédiction. Pour évaluer un modèle dans son ensemble, nous voulons mesurer sa performance moyenne sur toutes les données possibles, pas seulement sur les exemples que nous avons observés.

Pourquoi des variables aléatoires?¶

Une question naturelle se pose: si nous ajustons une fonction déterministe $f$ à des données, pourquoi avons-nous besoin de variables aléatoires et d’espérances? La fonction obtenue n’est-elle pas simplement une courbe fixe?

La réponse tient en un mot: généralisation. Nous ne nous intéressons pas vraiment à la performance sur les données d’entraînement car ces points sont déjà connus. Ce qui compte, c’est la performance sur des données futures que nous n’avons pas encore observées.

Considérons l’exemple de la distance de freinage. Les 50 mesures dans notre tableau sont un échantillon de toutes les mesures possibles. Si nous retournions sur le terrain et mesurions à nouveau, nous obtiendrions des valeurs légèrement différentes. En effet, le même véhicule à 20 mph ne s’arrête pas exactement à la même distance à chaque essai. Il y a de la variabilité intrinsèque: état de la route, température des freins, réflexes du conducteur.

Cette variabilité est capturée par une distribution de probabilité $p(x, y)$. Nos 50 points sont des tirages de cette distribution. La question fondamentale devient alors:

Notre modèle $f$ prédira-t-il bien sur de nouveaux tirages de cette même distribution?

La fonction $f$ elle-même est déterministe une fois entraînée. Mais son évaluation, savoir si elle prédit bien ou mal, dépend de quelles données futures elle rencontrera. Et ces données futures sont incertaines: elles seront tirées de $p(x, y)$, mais nous ne savons pas lesquelles.

Le risque formalise cette idée: c’est la perte moyenne que subira notre modèle $f$ lorsqu’il sera confronté à des données tirées de $p(x, y)$. C’est une mesure de performance prospective, tournée vers le futur.

Définition formelle¶

Le risque d’une fonction $f$ est l’espérance de la perte sur la distribution des données:

Décomposons cette formule étape par étape:

1. $\mathbb{E}_{(\mathbf{X},Y) \sim p}$: L’espérance mathématique signifie “moyenne sur tous les exemples possibles”. La notation $(\mathbf{X},Y) \sim p$ indique que nous tirons les paires $(\mathbf{x}, y)$ selon la distribution $p(\mathbf{x}, y)$ de la nature.

2. $\ell(Y, f(\mathbf{X}))$: Pour chaque exemple aléatoire $(\mathbf{X}, Y)$, nous calculons la perte entre la vraie valeur $Y$ et la prédiction $f(\mathbf{X})$ du modèle.

3. L’intégrale $\int \ell(y, f(\mathbf{x})) \, p(\mathbf{x}, y) \, d\mathbf{x} \, dy$: Cette intégrale calcule une moyenne pondérée. Pour chaque paire possible $(\mathbf{x}, y)$, nous multiplions la perte $\ell(y, f(\mathbf{x}))$ par la probabilité $p(\mathbf{x}, y)$ que cette paire apparaisse dans la nature, puis nous sommons (intégrons) sur toutes les paires possibles.

Exemple concret¶

Considérons un problème de classification binaire en 2D. Supposons que $\mathbf{x} \in [0, 1]^2$ et $y \in \{0, 1\}$. Pour calculer le risque, nous devrions:

1. Diviser l’espace $[0,1]^2$ en une grille fine (par exemple, $1000 \times 1000$ points)

2. Pour chaque point $\mathbf{x}$ de la grille, considérer les deux valeurs possibles de $y$ (0 et 1)

3. Pour chaque combinaison $(\mathbf{x}, y)$, calculer:

   • La probabilité $p(\mathbf{x}, y)$ que cette combinaison apparaisse

   • La perte $\ell(y, f(\mathbf{x}))$ si notre modèle prédit $f(\mathbf{x})$

4. Faire la somme pondérée: $\sum_{\mathbf{x}} \sum_{y \in \{0,1\}} \ell(y, f(\mathbf{x})) \cdot p(\mathbf{x}, y)$

En pratique, pour un espace continu, cette somme devient une intégrale sur un domaine continu, ce qui est encore plus complexe à calculer.

Visualisons ceci concrètement. La figure suivante montre un problème de classification binaire où chaque classe suit une distribution gaussienne en 2D. Les contours représentent la densité $p(x|y)$ pour chaque classe. La ligne pointillée est la frontière de décision d’un classificateur linéaire. Les régions ombrées indiquent où le classificateur fait des erreurs: la région rouge correspond aux points de classe 0 classés comme classe 1, et la région bleue correspond aux points de classe 1 classés comme classe 0.

Source

import numpy as np
import matplotlib.pyplot as plt

# Paramètres du mélange gaussien (classification binaire 2D)
mu0 = np.array([0.0, 0.0])
mu1 = np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])

def gaussian_pdf(x, mu, cov):
    """PDF d'une gaussienne multivariée."""
    d = len(mu)
    diff = x - mu
    cov_inv = np.linalg.inv(cov)
    mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
    return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))

# Create grid for visualization
x_range = np.linspace(-3, 5, 200)
y_range = np.linspace(-3, 4, 200)
X_grid, Y_grid = np.meshgrid(x_range, y_range)
pos = np.dstack([X_grid, Y_grid])

# Compute class-conditional densities
p_x_given_0 = gaussian_pdf(pos, mu0, cov)
p_x_given_1 = gaussian_pdf(pos, mu1, cov)

# Joint densities (with equal priors)
prior = 0.5
p_x_y0 = p_x_given_0 * (1 - prior)
p_x_y1 = p_x_given_1 * prior

# Linear decision boundary (Bayes optimal for equal covariances)
# w^T x + b = 0 where w = Sigma^{-1}(mu1 - mu0)
cov_inv = np.linalg.inv(cov)
w = cov_inv @ (mu1 - mu0)
b = -0.5 * (mu1 @ cov_inv @ mu1 - mu0 @ cov_inv @ mu0)

# Decision boundary: w[0]*x + w[1]*y + b = 0  =>  y = -(w[0]*x + b)/w[1]
x_boundary = np.linspace(-3, 5, 100)
y_boundary = -(w[0] * x_boundary + b) / w[1]

# Classifier prediction: classify as 1 if w^T x + b > 0
predictions = (w[0] * X_grid + w[1] * Y_grid + b) > 0

# Misclassification regions
# Class 0 misclassified as 1: true class is 0, but prediction is 1
misclass_0 = predictions  # region where we predict 1
# Class 1 misclassified as 0: true class is 1, but prediction is 0
misclass_1 = ~predictions  # region where we predict 0

fig, ax = plt.subplots(figsize=(8, 6))

# Plot class-conditional densities as contours
levels = [0.01, 0.05, 0.1, 0.15]
ax.contour(X_grid, Y_grid, p_x_given_0, levels=levels, colors='C0', alpha=0.7)
ax.contour(X_grid, Y_grid, p_x_given_1, levels=levels, colors='C1', alpha=0.7)

# Shade misclassification regions weighted by probability
# Red: class 0 points incorrectly classified as 1
error_region_0 = np.where(misclass_0, p_x_y0, 0)
# Blue: class 1 points incorrectly classified as 0  
error_region_1 = np.where(misclass_1, p_x_y1, 0)

ax.contourf(X_grid, Y_grid, error_region_0, levels=[0.001, 0.01, 0.05, 0.1], 
            colors=['#ff000010', '#ff000030', '#ff000050'], extend='max')
ax.contourf(X_grid, Y_grid, error_region_1, levels=[0.001, 0.01, 0.05, 0.1],
            colors=['#0000ff10', '#0000ff30', '#0000ff50'], extend='max')

# Decision boundary
ax.plot(x_boundary, y_boundary, 'k--', linewidth=2, label='Frontière de décision')

# Class centers
ax.scatter(*mu0, s=100, c='C0', marker='x', linewidths=3, zorder=5, label='Centre classe 0')
ax.scatter(*mu1, s=100, c='C1', marker='x', linewidths=3, zorder=5, label='Centre classe 1')

ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.legend(loc='upper left')
ax.set_aspect('equal')

plt.tight_layout()

Le risque est l’intégrale de la perte sur tout l’espace, pondérée par $p(\mathbf{x}, y)$. Les régions ombrées contribuent au risque: chaque point dans ces régions est mal classé, et sa contribution dépend de la densité de probabilité à cet endroit. Les régions denses proches de la frontière contribuent le plus au risque.

Pourquoi le risque est important¶

Le risque mesure ce que nous obtiendrons en moyenne si nous appliquons $f$ à de nouvelles données tirées de la même distribution. Un modèle avec un faible risque fait de bonnes prédictions en général, pas seulement sur les exemples d’entraînement. C’est exactement ce que nous voulons optimiser: un modèle qui performe bien sur des données jamais vues, pas seulement sur celles qu’il a déjà observées.

Cette quantité est ce que nous voulons minimiser. Le problème fondamental est que nous ne connaissons pas la distribution $p(\mathbf{x}, y)$ de la nature. Nous n’y avons accès qu’indirectement, via un échantillon fini $\mathcal{D}$.

Le risque empirique¶

Puisque le risque est inaccessible, nous l’approximons par une moyenne sur les données disponibles. Le risque empirique est:

Cette quantité est calculable: c’est la moyenne des pertes sur l’échantillon d’entraînement. Pour la perte 0-1, le risque empirique est le taux d’erreur sur les données d’entraînement. Pour la perte quadratique, c’est l’erreur quadratique moyenne.

Pourquoi le risque est-il inaccessible?¶

La nécessité d’utiliser le risque empirique découle de deux obstacles fondamentaux, l’un conceptuel et l’autre computationnel.

Obstacle 1: La distribution $p(\mathbf{x}, y)$ est inconnue¶

La nature possède une distribution $p(\mathbf{x}, y)$ qui génère les données, mais nous ne la connaissons pas. Nous n’observons qu’un échantillon fini $\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^N$ tiré de cette distribution.

L’ensemble $\mathcal{D}$ est une variable aléatoire: si nous répétions l’expérience de collecte de données, nous obtiendrions un échantillon différent. Cette perspective, adoptée notamment dans Murphy (2022), rappelle que nos conclusions dépendent de l’échantillon particulier que nous avons observé. C’est comme si nous regardions quelques gouttes d’eau d’un océan: nous pouvons analyser ces gouttes, mais un autre prélèvement donnerait des gouttes différentes.

Même si nous tentions d’estimer $p(\mathbf{x}, y)$ à partir des données (par exemple, via des techniques d’estimation de densité comme les mélanges de gaussiennes ou les estimateurs à noyau), nous n’obtiendrions qu’une approximation $\hat{p}(\mathbf{x}, y)$ de la vraie distribution. Cette approximation serait elle-même imparfaite et dépendrait de nos hypothèses sur la forme de la distribution.

La figure suivante illustre ce problème. À gauche, la vraie distribution $p(\mathbf{x}, y)$ que la nature utilise pour générer les données (que nous ne connaissons pas). À droite, un échantillon de $N = 50$ points tirés de cette distribution (ce que nous observons).

Source

import numpy as np
import matplotlib.pyplot as plt

# Paramètres du mélange gaussien
mu0, mu1 = np.array([0.0, 0.0]), np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])

def gaussian_pdf(x, mu, cov):
    d = len(mu)
    diff = x - mu
    cov_inv = np.linalg.inv(cov)
    mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
    return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))

# Générer un échantillon
rng = np.random.default_rng(42)
n = 50
X0 = rng.multivariate_normal(mu0, cov, n // 2)
X1 = rng.multivariate_normal(mu1, cov, n // 2)
X = np.vstack([X0, X1])
y = np.concatenate([np.zeros(n // 2), np.ones(n // 2)])

# Grille pour visualisation
x_range = np.linspace(-3, 5, 150)
y_range = np.linspace(-3, 4, 150)
X_grid, Y_grid = np.meshgrid(x_range, y_range)
pos = np.dstack([X_grid, Y_grid])

p_x_given_0 = gaussian_pdf(pos, mu0, cov)
p_x_given_1 = gaussian_pdf(pos, mu1, cov)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Left: True distribution (what nature knows)
ax = axes[0]
ax.contourf(X_grid, Y_grid, p_x_given_0, levels=15, cmap='Blues', alpha=0.6)
ax.contourf(X_grid, Y_grid, p_x_given_1, levels=15, cmap='Oranges', alpha=0.6)
ax.contour(X_grid, Y_grid, p_x_given_0, levels=5, colors='C0', alpha=0.8)
ax.contour(X_grid, Y_grid, p_x_given_1, levels=5, colors='C1', alpha=0.8)
ax.scatter(*mu0, s=100, c='C0', marker='x', linewidths=3, zorder=5)
ax.scatter(*mu1, s=100, c='C1', marker='x', linewidths=3, zorder=5)
ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('Distribution vraie $p(x, y)$\n(inconnue)')
ax.set_aspect('equal')

# Right: Finite sample (what we observe)
ax = axes[1]
ax.scatter(X[y == 0, 0], X[y == 0, 1], c='C0', alpha=0.7, s=50, label='Classe 0')
ax.scatter(X[y == 1, 0], X[y == 1, 1], c='C1', alpha=0.7, s=50, label='Classe 1')
ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title(f'Échantillon observé $\\mathcal{{D}}$\n($N = {len(X)}$ points)')
ax.legend()
ax.set_aspect('equal')

plt.tight_layout()

Nous ne voyons que les points à droite. La structure continue à gauche, incluant les contours, les densités, ainsi que les régions de haute et basse probabilité, nous est cachée. C’est à partir de ces quelques points que nous devons estimer la performance de notre modèle.

Obstacle 2: L’intégration est computationnellement intractable¶

Supposons, par un miracle, que nous connaissions exactement $p(\mathbf{x}, y)$. Pourrions-nous alors calculer le risque $\mathcal{R}(f) = \int \ell(y, f(\mathbf{x})) \, p(\mathbf{x}, y) \, d\mathbf{x} \, dy$?

La réponse est généralement non, pour plusieurs raisons:

Pour les espaces continus: L’intégrale est une intégrale de grande dimension. Si $\mathbf{x} \in \mathbb{R}^d$ avec $d$ grand (par exemple, $d = 1000$ pour des images ou $d = 10^6$ pour des données textuelles), nous devons intégrer sur un espace de dimension $d+1$.

Pour vous rappeler l’idée, en calcul on approche une intégrale en 1D par une somme: on découpe l’intervalle en petites tranches et on additionne des aires de rectangles ou de trapèzes. Par exemple, sur $[a,b]$:

Cette idée générale, qui consiste à remplacer une intégrale par une somme pondérée de valeurs de $g$ évaluées à des points $x_m$, s’appelle l’intégration numérique (ou quadrature).

Le problème en apprentissage est que notre intégrale n’est pas en 1D. Si on applique le même raisonnement en dimension $d$ en mettant, disons, $M$ points par dimension, on obtient une grille de taille $M^d$ (et ici $d$ peut être très grand). Le nombre de points à évaluer explose donc exponentiellement avec $d$. C’est exactement la malédiction de la dimensionnalité.

Pour les espaces discrets: Si $\mathbf{x}$ et $y$ sont discrets mais prennent de nombreuses valeurs, la somme $\sum_{\mathbf{x}} \sum_y \ell(y, f(\mathbf{x})) \cdot p(\mathbf{x}, y)$ peut avoir un nombre exponentiel de termes. Par exemple, si $\mathbf{x}$ est un vecteur binaire de dimension $d$, il y a $2^d$ valeurs possibles pour $\mathbf{x}$. Pour $d = 100$, cela fait déjà $2^{100} \approx 10^{30}$ termes à sommer, ce qui est computationnellement impossible.

Intégration de Monte Carlo: On pourrait penser utiliser l’intégration de Monte Carlo: tirer des échantillons $(\mathbf{x}, y)$ selon $p(\mathbf{x}, y)$ et estimer l’intégrale par la moyenne empirique. Mais pour obtenir une estimation précise du risque, nous aurions besoin d’un très grand nombre d’échantillons (potentiellement infini pour une précision parfaite). De plus, cela nécessiterait de pouvoir échantillonner efficacement depuis $p(\mathbf{x}, y)$, ce qui est lui-même un problème difficile si la distribution est complexe.

La figure suivante illustre la malédiction de la dimensionnalité. Avec seulement 10 points par dimension pour une quadrature numérique, le nombre total de points d’évaluation explose rapidement.

Source

import numpy as np
import matplotlib.pyplot as plt

# Number of grid points per dimension
points_per_dim = 10

# Dimensions to consider
dimensions = np.array([1, 2, 3, 5, 10, 20, 50, 100])

# Total grid points = points_per_dim^d
total_points = points_per_dim ** dimensions.astype(float)

fig, ax = plt.subplots(figsize=(8, 5))

bars = ax.bar(range(len(dimensions)), total_points, color='steelblue', edgecolor='black')

# Add reference lines
ax.axhline(y=1e9, color='C1', linestyle='--', alpha=0.7, label='1 milliard (limite pratique)')
ax.axhline(y=1e80, color='C3', linestyle=':', alpha=0.7, label='$10^{80}$ (atomes dans l\'univers)')

ax.set_yscale('log')
ax.set_xticks(range(len(dimensions)))
ax.set_xticklabels([f'd={d}' for d in dimensions])
ax.set_xlabel('Dimension de l\'espace des entrées')
ax.set_ylabel('Nombre de points de grille')
ax.set_title(f'Points nécessaires pour l\'intégration numérique\n({points_per_dim} points par dimension)')
ax.legend(loc='upper left')

# Annotate a few bars
for i, (d, n) in enumerate(zip(dimensions, total_points)):
    if d <= 5:
        ax.annotate(f'$10^{{{d}}}$', (i, n), ha='center', va='bottom', fontsize=9)
    elif d == 10:
        ax.annotate(f'$10^{{{10}}}$', (i, n), ha='center', va='bottom', fontsize=9)
    elif d == 100:
        ax.annotate(f'$10^{{{100}}}$', (i, n), ha='center', va='bottom', fontsize=9)

ax.set_ylim(1, 1e105)

plt.tight_layout()

En dimension 10, il faut déjà 10¹⁰ points, soit dix milliards. En dimension 100, il en faut 10¹⁰⁰, un nombre qui dépasse le nombre d’atomes dans l’univers observable. L’intégration numérique directe est donc impossible pour les problèmes de haute dimension, même si nous connaissions $p(\mathbf{x}, y)$ exactement.

Le risque empirique comme seule option pratique¶

Face à ces obstacles, le risque empirique est notre seule option calculable. Mais il y a une bonne nouvelle: le risque empirique est une forme d’intégration de Monte Carlo, et Monte Carlo a une propriété remarquable.

La complexité de Monte Carlo est indépendante de la dimension $d$. Elle ne dépend que du nombre d’échantillons $N$. C’est cette propriété qui rend l’apprentissage possible en haute dimension. De plus, nous n’avons pas besoin de connaître la valeur numérique de $p(\mathbf{x},y)$: nous avons seulement besoin de pouvoir tirer des échantillons de cette distribution. C’est exactement ce que nos données d’entraînement nous fournissent.

Le risque empirique remplace l’intégrale sur la distribution inconnue par une moyenne sur l’échantillon fini que nous possédons:

Cette formule est directe à évaluer: nous parcourons nos $N$ exemples d’entraînement, calculons la perte pour chacun, et faisons la moyenne.

Reprenons les données de freinage. Divisons-les en deux parties: les mesures à vitesses faibles (4-19 mph) pour l’entraînement, et les mesures à vitesses élevées (20-25 mph) pour le test. Le risque empirique sur l’ensemble d’entraînement mesure la qualité de l’ajustement. Le risque empirique sur l’ensemble de test estime la performance sur des vitesses non vues pendant l’entraînement.

Source

import numpy as np
import matplotlib.pyplot as plt

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

# Split: train on low speeds, test on high speeds
train_mask = speed < 20
test_mask = speed >= 20

speed_train, dist_train = speed[train_mask], dist[train_mask]
speed_test, dist_test = speed[test_mask], dist[test_mask]

# Fit on training data
coeffs = np.polyfit(speed_train, dist_train, 2)

# Compute MSE on train and test
pred_train = np.polyval(coeffs, speed_train)
pred_test = np.polyval(coeffs, speed_test)
mse_train = np.mean((dist_train - pred_train)**2)
mse_test = np.mean((dist_test - pred_test)**2)

plt.figure(figsize=(7, 4))
plt.scatter(speed_train, dist_train, alpha=0.7, label=f'Entraînement (EQM={mse_train:.1f})')
plt.scatter(speed_test, dist_test, alpha=0.7, marker='s', label=f'Test (EQM={mse_test:.1f})')

speed_grid = np.linspace(4, 28, 100)
plt.plot(speed_grid, np.polyval(coeffs, speed_grid), 'k--', alpha=0.6, label='Fonction ajustée')

plt.axvline(x=20, color='gray', linestyle=':', alpha=0.5)
plt.xlabel('Vitesse (mph)')
plt.ylabel('Distance (ft)')
plt.legend()
plt.tight_layout()

Dans cet exemple, l’EQM (MSE, erreur quadratique moyenne) sur l’ensemble de test est plus élevé que sur l’ensemble d’entraînement. Cet écart est typique: la fonction a été optimisée pour les données d’entraînement, pas pour les données de test.

Sous l’hypothèse que les exemples $(\mathbf{x}_i, y_i)$ sont tirés indépendamment et identiquement distribués (i.i.d.) selon $p(\mathbf{x}, y)$, le risque empirique est un estimateur non biaisé du vrai risque: $\mathbb{E}[\hat{\mathcal{R}}(f, \mathcal{D})] = \mathcal{R}(f)$. Cela signifie qu’en moyenne, sur tous les échantillons possibles, le risque empirique est égal au vrai risque.

Par la loi des grands nombres, lorsque $N \to \infty$, le risque empirique converge vers le vrai risque (presque sûrement). Avec suffisamment de données, si l’échantillon est représentatif de la distribution, le risque empirique devrait être proche du risque.

La figure suivante illustre cette convergence. Nous utilisons le problème de classification gaussienne pour lequel nous pouvons calculer le vrai risque analytiquement. Chaque courbe montre l’évolution du risque empirique pour un échantillon de taille croissante. Toutes les courbes convergent vers le vrai risque (ligne pointillée), mais avec des fluctuations qui diminuent à mesure que $N$ augmente.

Source

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Paramètres du mélange gaussien
mu0, mu1 = np.array([0.0, 0.0]), np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])

def gaussian_pdf(x, mu, cov):
    d = len(mu)
    diff = x - mu
    cov_inv = np.linalg.inv(cov)
    mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
    return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))

# Compute Bayes-optimal classifier error rate (true risk)
# For Gaussian classes with equal covariance, the Bayes error is:
# P(error) = Phi(-d/2) where d is the Mahalanobis distance between means
cov_inv = np.linalg.inv(cov)
d_squared = (mu1 - mu0) @ cov_inv @ (mu1 - mu0)
d = np.sqrt(d_squared)
true_risk = norm.cdf(-d / 2)

# Simulate empirical risk for different sample sizes
sample_sizes = np.arange(10, 1001, 10)
n_runs = 20

fig, ax = plt.subplots(figsize=(9, 5))

# Store all runs for confidence band
all_risks = np.zeros((n_runs, len(sample_sizes)))

for run in range(n_runs):
    empirical_risks = []
    # Generate a large dataset and compute cumulative empirical risk
    rng = np.random.default_rng(run)
    n_total = 1000
    X0 = rng.multivariate_normal(mu0, cov, n_total // 2)
    X1 = rng.multivariate_normal(mu1, cov, n_total // 2)
    X = np.vstack([X0, X1])
    y = np.concatenate([np.zeros(n_total // 2), np.ones(n_total // 2)])
    perm = rng.permutation(n_total)
    X, y = X[perm], y[perm]
    
    # Bayes-optimal classifier: predict 1 if w^T x + b > 0
    w = cov_inv @ (mu1 - mu0)
    b = -0.5 * (mu1 @ cov_inv @ mu1 - mu0 @ cov_inv @ mu0)
    
    for n in sample_sizes:
        X_n, y_n = X[:n], y[:n]
        predictions = (X_n @ w + b > 0).astype(float)
        emp_risk = np.mean(predictions != y_n)
        empirical_risks.append(emp_risk)
    
    all_risks[run] = empirical_risks
    ax.plot(sample_sizes, empirical_risks, 'C0-', alpha=0.15, linewidth=0.8)

# Mean and confidence bands
mean_risk = np.mean(all_risks, axis=0)
std_risk = np.std(all_risks, axis=0)
ax.fill_between(sample_sizes, mean_risk - 2*std_risk, mean_risk + 2*std_risk, 
                alpha=0.3, color='C0', label='Intervalle ± 2 écarts-types')
ax.plot(sample_sizes, mean_risk, 'C0-', linewidth=2, label='Moyenne empirique')

# True risk
ax.axhline(y=true_risk, color='C3', linestyle='--', linewidth=2, 
           label=f'Vrai risque = {true_risk:.3f}')

ax.set_xlabel('Taille de l\'échantillon $N$')
ax.set_ylabel('Risque empirique (taux d\'erreur)')
ax.set_title('Convergence du risque empirique vers le vrai risque')
ax.legend(loc='upper right')
ax.set_xlim(0, 1000)
ax.set_ylim(0, 0.35)

plt.tight_layout()

Avec $N = 50$, le risque empirique peut facilement varier de 0.10 à 0.25 selon l’échantillon. Avec $N = 500$, la variabilité est beaucoup plus faible. C’est la loi des grands nombres en action: plus l’échantillon est grand, plus l’estimation est précise.

Le compromis fondamental¶

Cette situation crée un compromis fondamental en apprentissage automatique:

• Ce que nous voulons minimiser: Le risque $\mathcal{R}(f)$, qui mesure la performance sur toutes les données possibles

• Ce que nous pouvons minimiser: Le risque empirique $\hat{\mathcal{R}}(f, \mathcal{D})$, qui mesure la performance sur nos données d’entraînement

L’écart entre ces deux quantités est au cœur de l’apprentissage automatique. Un modèle peut avoir un risque empirique très faible (il performe bien sur les données d’entraînement) tout en ayant un risque élevé (il performe mal sur de nouvelles données). C’est le problème du surapprentissage, que nous explorerons plus en détail dans le chapitre sur la généralisation.

La question de savoir quand et à quelle vitesse l’approximation du risque par le risque empirique est fiable relève de la théorie de la généralisation, que nous aborderons au chapitre suivant.

Minimisation du risque empirique¶

Nous avons maintenant les éléments pour formuler l’apprentissage comme un problème d’optimisation. Nous cherchons la fonction $f$ dans une classe $\mathcal{F}$ qui minimise le risque:

Puisque le risque est inaccessible, nous le remplaçons par le risque empirique:

Ce principe est la minimisation du risque empirique (MRE): choisir la fonction qui fait le moins d’erreurs sur les données d’entraînement, en espérant que cette performance se transfère aux nouvelles données.

La classe $\mathcal{F}$ est notre classe d’hypothèses. Elle représente l’ensemble des fonctions que nous sommes prêts à considérer. Le choix de $\mathcal{F}$ encode nos hypothèses sur la forme de la relation entre entrées et sorties.

Un premier exemple: les modèles linéaires¶

Pour rendre ces concepts concrets, considérons la classe la plus simple: les modèles linéaires. Un modèle linéaire suppose que la sortie est une combinaison linéaire des entrées:

où $\mathbf{x} \in \mathbb{R}^{d+1}$ est le vecteur d’entrée augmenté d’un 1 pour le biais ($x_0 = 1$), et $\boldsymbol{\theta} \in \mathbb{R}^{d+1}$ est le vecteur de paramètres contenant le biais $\theta_0$ et les poids $\theta_1, \ldots, \theta_d$.

Cette forme est restrictive: elle suppose que la relation entre entrées et sorties est linéaire. Pour les données de freinage, cela signifierait que la distance est proportionnelle à la vitesse, ce qui n’est pas le cas (la relation est plutôt quadratique). Néanmoins, les modèles linéaires sont utiles comme point de départ, et nous verrons comment les étendre pour capturer des relations non linéaires.

Avec cette classe $\mathcal{F}$ fixée, l’apprentissage consiste à trouver les paramètres $\boldsymbol{\theta}$ qui minimisent le risque empirique. Pour la perte quadratique, cela revient à minimiser la somme des carrés des résidus.

Mais quand le minimiseur du risque empirique a-t-il un faible risque? Si $\hat{f}$ minimise $\hat{\mathcal{R}}$ et $f^\star$ minimise $\mathcal{R}$, nous voulons que $\mathcal{R}(\hat{f})$ soit proche de $\mathcal{R}(f^\star)$. La réponse dépend de la taille de l’échantillon $N$, de la complexité de la classe $\mathcal{F}$, et de propriétés de la distribution $p$.

Résoudre le problème d’optimisation¶

Quand nous écrivons $\arg\min_{\boldsymbol{\theta}}$, nous cherchons les paramètres qui rendent la fonction objectif aussi petite que possible. Mais comment trouver ces paramètres en pratique?

La réponse dépend de la forme du problème:

• Pour certains problèmes, comme la régression linéaire avec perte quadratique, nous pouvons dériver une solution analytique en posant le gradient égal à zéro et en résolvant le système d’équations résultant.

• Pour d’autres, nous devons recourir à des algorithmes itératifs (comme la descente de gradient) ou à des solveurs spécialisés (comme la programmation quadratique pour les SVM).

Exemple: solution analytique pour la régression linéaire (MCO)¶

Pour illustrer les solutions analytiques, considérons la régression linéaire avec perte quadratique. L’objectif est de minimiser la somme des carrés des résidus:

où $\mathbf{X}$ est la matrice $N \times (d+1)$ des entrées (avec une colonne de 1 pour le biais) et $\mathbf{y}$ est le vecteur des sorties.

En développant et en calculant le gradient:

En posant le gradient égal à zéro, nous obtenons les équations normales:

Si la matrice $\mathbf{X}^\top \mathbf{X}$ est inversible, la solution unique est:

Cette solution porte le nom d’estimateur des moindres carrés ordinaires (MCO, ou ordinary least squares, OLS). Elle peut être calculée directement sans itération, ce qui en fait un exemple classique de solution analytique.

Le chapitre sur l’optimisation présentera ces méthodes en détail. Pour l’instant, nous nous concentrons sur la formulation des problèmes d’apprentissage comme problèmes d’optimisation, en gardant à l’esprit que des outils existent pour les résoudre.

Généralisation¶

La différence entre le risque et le risque empirique est l’écart de généralisation:

Un modèle qui minimise le risque empirique peut avoir un risque élevé si cet écart est grand. Ce phénomène est le surapprentissage: le modèle s’ajuste aux particularités de l’échantillon d’entraînement, y compris le bruit, plutôt qu’aux régularités sous-jacentes. L’erreur d’entraînement est faible, mais l’erreur sur de nouvelles données est élevée.

À l’inverse, un modèle trop simple peut avoir un risque empirique et un risque tous deux élevés. C’est le sous-apprentissage: le modèle n’a pas la capacité de capturer la structure des données.

Extrapolation¶

Un cas particulier de mauvaise généralisation est l’extrapolation: prédire pour des entrées en dehors de la plage des données d’entraînement. Même un modèle bien ajusté peut échouer spectaculairement lorsqu’on lui demande de prédire au-delà de ce qu’il a vu.

Considérons des essais en soufflerie pour mesurer la portance d’une aile à différentes vitesses. Les tests sont effectués entre 20 et 60 m/s. L’ingénieur veut prédire la portance à 100 m/s.

Source

import numpy as np
import matplotlib.pyplot as plt

# Données de portance aérodynamique (simulées)
np.random.seed(42)
rho, S, C_L = 1.225, 20.0, 0.5
v_train = np.linspace(20, 60, 8)
L_true_train = 0.5 * rho * v_train**2 * S * C_L
L_train = L_true_train + np.random.normal(0, 400, len(v_train))

coeffs_2 = np.polyfit(v_train, L_train, 2)
coeffs_5 = np.polyfit(v_train, L_train, 5)

v_extrap = np.linspace(15, 110, 200)
L_true_extrap = 0.5 * rho * v_extrap**2 * S * C_L

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

for ax, coeffs, deg in zip(axes, [coeffs_2, coeffs_5], [2, 5]):
    ax.scatter(v_train, L_train, s=50, zorder=5, label='Observations')
    ax.plot(v_extrap, L_true_extrap, 'g-', alpha=0.5, label='Vraie relation')
    L_pred = np.polyval(coeffs, v_extrap)
    ax.plot(v_extrap, L_pred, 'k--', label=f'Polynôme degré {deg}')
    ax.axvline(60, color='gray', linestyle=':', alpha=0.5)
    ax.axvspan(60, 110, alpha=0.1, color='red')
    ax.set_xlabel('Vitesse (m/s)')
    ax.set_ylabel('Portance (N)')
    ax.set_title(f'Degré {deg}')
    ax.legend(loc='upper left')
    ax.set_ylim(-5000, 80000)
    ax.text(85, 5000, 'Extrapolation', ha='center', fontsize=10, color='red', alpha=0.7)

plt.tight_layout()

Le polynôme de degré 2 (qui correspond au vrai modèle physique $L \propto v^2$) extrapole correctement. Le polynôme de degré 5, bien qu’il ajuste aussi bien les données d’entraînement, diverge complètement en dehors de la plage observée.

Régularisation¶

Une manière de contrôler le surapprentissage consiste à pénaliser la complexité du modèle directement dans la fonction objectif. Le risque empirique régularisé est:

où $C(\boldsymbol{\theta})$ mesure la complexité du modèle et $\lambda \geq 0$ contrôle l’intensité de la pénalisation. Un choix courant est la régularisation $\ell_2$ (ou weight decay):

Cette pénalisation pousse les paramètres vers zéro, ce qui a pour effet de lisser la fonction apprise. En régression linéaire, l’ajout de cette pénalité donne la régression ridge:

Illustrons l’effet de la régularisation sur le même problème de régression polynomiale. Avec un polynôme de degré 15 et différentes valeurs de $\lambda$:

Source

import numpy as np
import matplotlib.pyplot as plt

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

# Train/test split
np.random.seed(42)
indices = np.random.permutation(len(speed))
train_idx, test_idx = indices[:35], indices[35:]
speed_train, dist_train = speed[train_idx], dist[train_idx]
speed_test, dist_test = speed[test_idx], dist[test_idx]

# Build polynomial features (degree 15)
degree = 15
def poly_features(x, deg):
    return np.vstack([x**i for i in range(deg+1)]).T

X_train = poly_features(speed_train, degree)
X_test = poly_features(speed_test, degree)

# Ridge regression for different lambda values
lambdas = [0, 1e-6, 1e-3, 1]
fig, axes = plt.subplots(2, 2, figsize=(10, 8))

for ax, lam in zip(axes.flat, lambdas):
    # Ridge solution: (X^T X + lambda I)^{-1} X^T y
    I = np.eye(X_train.shape[1])
    I[0, 0] = 0  # Don't regularize bias
    w = np.linalg.solve(X_train.T @ X_train + lam * I, X_train.T @ dist_train)
    
    # Predictions
    pred_train = X_train @ w
    pred_test = X_test @ w
    mse_train = np.mean((dist_train - pred_train)**2)
    mse_test = np.mean((dist_test - pred_test)**2)
    
    # Plot
    ax.scatter(speed_train, dist_train, alpha=0.6, s=30, label='Entraînement')
    ax.scatter(speed_test, dist_test, alpha=0.6, s=30, marker='s', label='Test')
    
    speed_grid = np.linspace(3, 26, 200)
    X_grid = poly_features(speed_grid, degree)
    pred_grid = X_grid @ w
    pred_grid = np.clip(pred_grid, -50, 200)
    ax.plot(speed_grid, pred_grid, 'k-', alpha=0.7)
    
    ax.set_xlim(3, 26)
    ax.set_ylim(-20, 150)
    ax.set_xlabel('Vitesse (mph)')
    ax.set_ylabel('Distance (ft)')
    ax.set_title(f'$\\lambda$ = {lam}: Entr. EQM={mse_train:.1f}, Test EQM={mse_test:.1f}')
    if lam == 0:
        ax.legend()

plt.tight_layout()

Sans régularisation ($\lambda = 0$), le polynôme de degré 15 oscille fortement. Avec une régularisation modérée ($\lambda = 10^{-3}$), les oscillations sont atténuées et l’erreur de test diminue. Avec une régularisation trop forte ($\lambda = 1$), le modèle devient trop contraint et sous-apprend.

Solution analytique de la régression ridge¶

Comme pour les moindres carrés ordinaires, la régression ridge admet une solution analytique. L’objectif régularisé est:

En développant et en calculant le gradient:

En posant le gradient égal à zéro, nous obtenons les équations normales régularisées:

La solution est:

Comparons avec la solution MCO: $\hat{\boldsymbol{\theta}}_{\text{MCO}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}$. La seule différence est l’ajout du terme $\lambda \mathbf{I}$ à la matrice $\mathbf{X}^\top \mathbf{X}$.

Solution via décomposition en valeurs singulières (SVD)¶

Cette section présente une autre façon d’exprimer les solutions MCO et Ridge, en utilisant la décomposition en valeurs singulières (SVD). Si vous n’avez jamais rencontré la SVD, ne vous inquiétez pas: nous allons l’introduire progressivement. Cette approche n’est pas strictement nécessaire pour comprendre la régression ridge, mais elle offre une interprétation géométrique très éclairante qui révèle pourquoi la régularisation fonctionne.

Qu’est-ce que la SVD? Si vous avez déjà rencontré la décomposition en valeurs propres, la SVD en est une généralisation. Pour une matrice carrée symétrique $\mathbf{A}$, la décomposition en valeurs propres s’écrit $\mathbf{A} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^\top$, où $\mathbf{Q}$ contient les vecteurs propres et $\boldsymbol{\Lambda}$ les valeurs propres. La SVD généralise cette idée à n’importe quelle matrice, même rectangulaire.

Pour une matrice $\mathbf{X}$ de données, la SVD la réécrit comme le produit de trois matrices:

Lien avec la décomposition en valeurs propres: Les colonnes de $\mathbf{V}$ sont les vecteurs propres de $\mathbf{X}^\top \mathbf{X}$, et les valeurs singulières $d_j$ sont les racines carrées des valeurs propres de $\mathbf{X}^\top \mathbf{X}$. Autrement dit, si $\mathbf{X}^\top \mathbf{X} = \mathbf{V} \boldsymbol{\Lambda} \mathbf{V}^\top$ est la décomposition en valeurs propres de $\mathbf{X}^\top \mathbf{X}$, alors $d_j = \sqrt{\lambda_j}$ où $\lambda_j$ sont les valeurs propres. De même, les colonnes de $\mathbf{U}$ sont les vecteurs propres de $\mathbf{X} \mathbf{X}^\top$.

Cette connexion est utile car $\mathbf{X}^\top \mathbf{X}$ apparaît naturellement dans la régression (c’est la matrice que nous inversons pour MCO). Les valeurs singulières $d_j$ nous renseignent donc directement sur le “conditionnement” de cette matrice: si certaines valeurs singulières sont très petites, alors $\mathbf{X}^\top \mathbf{X}$ est proche d’être singulière (non inversible).

Interprétation géométrique:

• $\mathbf{V}$ contient les directions principales dans l’espace des caractéristiques (les colonnes $\mathbf{v}_j$ sont orthonormales). Ces directions correspondent aux axes le long desquels la matrice $\mathbf{X}$ transforme les vecteurs de manière la plus efficace.

• $\mathbf{D}$ est une matrice diagonale contenant les valeurs singulières $d_1 \geq d_2 \geq \cdots \geq d_d \geq 0$, ordonnées du plus grand au plus petit. Chaque valeur singulière $d_j$ mesure l’amplitude de la transformation le long de la direction $\mathbf{v}_j$. Une grande valeur singulière signifie que la transformation est forte dans cette direction; une petite valeur singulière signifie que la transformation est faible.

• $\mathbf{U}$ contient les directions correspondantes dans l’espace des observations (les colonnes $\mathbf{u}_j$ sont orthonormales). Chaque $\mathbf{u}_j$ indique comment les observations se projettent sur la direction principale $\mathbf{v}_j$.

Solution MCO via SVD: En utilisant cette décomposition, la solution MCO peut s’écrire:

Cette formule décompose la solution en une somme de contributions le long de chaque direction principale $\mathbf{v}_j$. Le terme $\frac{\mathbf{u}_j^\top \mathbf{y}}{d_j}$ mesure combien la sortie $\mathbf{y}$ s’aligne avec la direction $\mathbf{u}_j$, divisé par l’amplitude $d_j$ de cette direction. Notez que diviser par une petite valeur singulière $d_j$ peut amplifier le bruit, ce qui explique pourquoi MCO peut être instable quand certaines directions ont de petites valeurs singulières.

Solution Ridge via SVD: Pour Ridge, la solution devient:

La différence avec MCO est le facteur de rétrécissement $\frac{d_j^2}{d_j^2 + \lambda}$ qui multiplie chaque terme. Ce facteur est toujours inférieur à 1, ce qui “rétrécit” chaque composante vers zéro. L’effet clé est que ce rétrécissement est différencié: les directions avec de petites valeurs singulières sont rétrécies plus fortement que celles avec de grandes valeurs singulières.

Interprétation géométrique: Les directions principales $\mathbf{v}_j$ définissent les axes d’une ellipse de confiance dans l’espace des coefficients. Pour MCO, les longueurs des demi-axes sont proportionnelles à $1/d_j$: plus une direction a une petite valeur singulière $d_j$, plus la variance d’estimation est grande (nous verrons cette distinction cruciale plus loin). Avec Ridge, ces longueurs deviennent proportionnelles à $\frac{d_j}{d_j^2 + \lambda}$: l’ellipse se rétrécit globalement, et plus rapidement le long des directions associées aux petites valeurs singulières. C’est exactement ce que nous voulons: réduire la variance d’estimation là où elle est la plus grande, c’est-à-dire dans les directions où les données sont peu dispersées.

Avantages numériques: Au-delà de l’interprétation, la SVD offre aussi des avantages pratiques. Elle est plus stable numériquement que l’inversion directe de $\mathbf{X}^\top \mathbf{X}$, surtout quand cette matrice est mal conditionnée (c’est-à-dire quand certaines valeurs singulières sont très petites). Les algorithmes SVD gèrent mieux ces cas délicats.

Visualisation: ellipse des données et vecteurs singuliers¶

Pour rendre ces concepts concrets, visualisons ce que la SVD capture sur un nuage de données 2D. Générons des points suivant une distribution gaussienne avec une covariance non triviale (les deux variables sont corrélées).

Source

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

np.random.seed(42)

# Générer des données gaussiennes corrélées
n_points = 200
mean = [0, 0]
cov = [[2.0, 1.2], [1.2, 1.0]]  # Covariance non diagonale
X = np.random.multivariate_normal(mean, cov, n_points)

# Centrer les données
X_centered = X - X.mean(axis=0)

# SVD de la matrice de données centrée
U, d, Vt = np.linalg.svd(X_centered, full_matrices=False)
V = Vt.T

# Les valeurs singulières sont liées aux écarts-types: d_j / sqrt(N-1)
# Pour l'ellipse, nous utilisons les écarts-types dans chaque direction
std_1 = d[0] / np.sqrt(n_points - 1)
std_2 = d[1] / np.sqrt(n_points - 1)

# Créer la figure
fig, ax = plt.subplots(figsize=(8, 6))

# Tracer les points
ax.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.5, s=20, c='tab:blue', label='Données')

# Tracer les vecteurs singuliers (directions principales)
origin = [0, 0]
scale = 2  # Facteur d'échelle pour la visualisation

# Premier vecteur singulier (direction de plus grande variance)
ax.annotate('', xy=V[:, 0] * std_1 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:red', lw=2.5))
ax.annotate('', xy=-V[:, 0] * std_1 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:red', lw=2.5))

# Deuxième vecteur singulier (direction de plus petite variance)
ax.annotate('', xy=V[:, 1] * std_2 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:orange', lw=2.5))
ax.annotate('', xy=-V[:, 1] * std_2 * scale, xytext=origin,
            arrowprops=dict(arrowstyle='->', color='tab:orange', lw=2.5))

# Ellipse de confiance (2 écarts-types)
angle = np.degrees(np.arctan2(V[1, 0], V[0, 0]))
ellipse = Ellipse(xy=(0, 0), width=4*std_1, height=4*std_2, angle=angle,
                  fill=False, edgecolor='gray', linestyle='--', linewidth=1.5)
ax.add_patch(ellipse)

# Annotations
ax.text(V[0, 0] * std_1 * scale * 1.15, V[1, 0] * std_1 * scale * 1.15, 
        f'$\\mathbf{{v}}_1$ ($d_1 = {d[0]:.1f}$)', fontsize=11, color='tab:red')
ax.text(V[0, 1] * std_2 * scale * 1.3, V[1, 1] * std_2 * scale * 1.3, 
        f'$\\mathbf{{v}}_2$ ($d_2 = {d[1]:.1f}$)', fontsize=11, color='tab:orange')

ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('Nuage gaussien et directions principales (SVD)')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.5)
ax.axvline(0, color='gray', linewidth=0.5)
ax.set_xlim(-4, 4)
ax.set_ylim(-3, 3)

plt.tight_layout()

La figure montre un nuage de 200 points tirés d’une gaussienne 2D. Les flèches représentent les vecteurs singuliers $\mathbf{v}_1$ et $\mathbf{v}_2$:

• Le vecteur $\mathbf{v}_1$ (rouge) pointe dans la direction de plus grande variance. La valeur singulière $d_1$ mesure l’amplitude de la dispersion dans cette direction.

• Le vecteur $\mathbf{v}_2$ (orange) pointe dans la direction de plus petite variance, perpendiculaire à $\mathbf{v}_1$. La valeur singulière $d_2$ est plus petite.

L’ellipse en pointillés représente la région contenant environ 95% des données si elles suivent exactement la distribution gaussienne. Ses axes coïncident avec les vecteurs singuliers, et les longueurs des demi-axes sont proportionnelles aux valeurs singulières.

Cette visualisation illustre pourquoi la SVD est si utile: elle identifie automatiquement les axes naturels des données. Dans le contexte de la régression, si les caractéristiques forment un nuage allongé (valeurs singulières très différentes), alors certaines directions contiennent beaucoup d’information (grandes valeurs singulières) tandis que d’autres en contiennent peu (petites valeurs singulières).

Deux variances: données vs estimation¶

Avant d’aller plus loin, clarifions une source fréquente de confusion. Le mot variance apparaît dans deux contextes très différents lorsqu’on parle de SVD et de régularisation:

1. Variance des données (dispersion): mesure l’étalement des données le long d’une direction $\mathbf{v}_j$. Elle est proportionnelle à $d_j^2$. Une grande valeur singulière $d_j$ signifie que les données sont très dispersées dans cette direction.

2. Variance d’estimation (incertitude): mesure l’incertitude sur notre estimé $\hat{\theta}_j$ du paramètre correspondant à la direction $\mathbf{v}_j$. Elle est proportionnelle à $1/d_j^2$. Une petite valeur singulière $d_j$ signifie que notre estimé est très incertain.

Ces deux variances sont inversement reliées:

Intuition: Imaginez estimer une pente à partir de données. Si les points sont très étalés horizontalement (grande variance des données en $x$), la pente est facile à déterminer avec précision (faible variance d’estimation). Si les points sont tous regroupés (petite variance des données), la pente est très incertaine (grande variance d’estimation).

C’est cette relation inverse qui explique le comportement de Ridge:

• Ridge rétrécit les directions où $d_j$ est petit (faible variance des données)

• Ce sont précisément les directions où la variance d’estimation est grande

• En rétrécissant ces directions, Ridge réduit la variance d’estimation au prix d’un biais

Ainsi, quand nous disons que « Ridge contrôle la variance », nous parlons de la variance d’estimation des paramètres, pas de la variance des données. La régularisation n’affecte pas la dispersion des données; elle réduit l’incertitude de nos estimés en les « tirant » vers zéro.

Spectre des valeurs singulières et rang effectif¶

En pratique, les données réelles ont souvent des dizaines ou des centaines de dimensions. Comment se comportent les valeurs singulières dans ce cas? Examinons un exemple avec des données de dimension plus élevée.

Source

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

# Simuler des données avec une structure de rang bas + bruit
n_samples = 100
n_features = 30

# Vraie structure: combinaison de 5 facteurs latents
n_latent = 5
latent = np.random.randn(n_samples, n_latent)
loadings = np.random.randn(n_latent, n_features)
X_signal = latent @ loadings

# Ajouter du bruit
noise_level = 0.5
X = X_signal + noise_level * np.random.randn(n_samples, n_features)

# Centrer
X_centered = X - X.mean(axis=0)

# SVD
U, d, Vt = np.linalg.svd(X_centered, full_matrices=False)

# Variance expliquée
variance_explained = d**2 / np.sum(d**2)
cumulative_variance = np.cumsum(variance_explained)

# Créer la figure
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# Panneau gauche: spectre des valeurs singulières
ax1 = axes[0]
ax1.semilogy(range(1, len(d)+1), d, 'o-', markersize=6, linewidth=1.5, color='tab:blue')
ax1.axhline(d[n_latent], color='tab:red', linestyle='--', linewidth=1.5, 
            label=f'Seuil (rang effectif = {n_latent})')
ax1.fill_between(range(1, n_latent+1), d[:n_latent], alpha=0.3, color='tab:green', label='Signal')
ax1.fill_between(range(n_latent+1, len(d)+1), d[n_latent:], alpha=0.3, color='tab:orange', label='Bruit')
ax1.set_xlabel('Indice $j$')
ax1.set_ylabel('Valeur singulière $d_j$ (échelle log)')
ax1.set_title('Spectre des valeurs singulières')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim(0.5, len(d)+0.5)

# Panneau droit: variance expliquée cumulative
ax2 = axes[1]
ax2.plot(range(1, len(d)+1), cumulative_variance * 100, 'o-', markersize=6, 
         linewidth=1.5, color='tab:blue')
ax2.axhline(95, color='gray', linestyle='--', linewidth=1, label='Seuil 95%')
ax2.axvline(n_latent, color='tab:red', linestyle='--', linewidth=1.5)

# Trouver k pour 95% de variance
k_95 = np.searchsorted(cumulative_variance, 0.95) + 1
ax2.scatter([k_95], [cumulative_variance[k_95-1]*100], s=100, color='tab:red', zorder=5)
ax2.annotate(f'{k_95} composantes\npour 95%', xy=(k_95, cumulative_variance[k_95-1]*100),
             xytext=(k_95+5, cumulative_variance[k_95-1]*100-10), fontsize=10,
             arrowprops=dict(arrowstyle='->', color='gray'))

ax2.set_xlabel('Nombre de composantes $k$')
ax2.set_ylabel('Variance expliquée cumulative (%)')
ax2.set_title('Variance expliquée')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0.5, len(d)+0.5)
ax2.set_ylim(0, 105)

plt.tight_layout()

Le panneau de gauche montre le spectre des valeurs singulières en échelle logarithmique. On observe un schéma typique:

• Les premières valeurs singulières sont grandes: elles correspondent aux directions du signal, la vraie structure sous-jacente des données.

• Après un certain point (ici, autour de $j = 5$), les valeurs singulières chutent et forment un “plancher”: ce sont les directions du bruit.

Le rang effectif est le nombre de valeurs singulières significativement au-dessus du plancher de bruit. Dans cet exemple, nous avons simulé 5 facteurs latents, et le spectre révèle bien cette structure: les 5 premières valeurs singulières dominent.

Le panneau de droite montre la variance expliquée cumulative. C’est un outil pratique pour choisir combien de composantes retenir:

• Critère du seuil: Retenir assez de composantes pour expliquer 95% (ou 99%) de la variance.

• Critère du coude: Chercher le “coude” dans le spectre où les valeurs singulières cessent de décroître rapidement.

• Critère du gap: Si le spectre présente un saut net (comme ici entre $d_5$ et $d_6$), c’est un bon point de coupure.

De la troncature à la réduction de dimension (ACP)¶

L’analyse ci-dessus suggère une idée: si les dernières directions ne contiennent que du bruit, pourquoi ne pas simplement les ignorer? Au lieu de rétrécir les coefficients comme Ridge, nous pourrions tronquer la représentation en ne gardant que les $k$ premières directions principales.

C’est exactement l’idée de l’analyse en composantes principales (ACP). Au lieu de travailler avec les $d$ caractéristiques originales, nous projetons les données sur les $k$ premiers vecteurs singuliers:

où $\mathbf{V}_k$ contient les $k$ premiers vecteurs singuliers (les colonnes de $\mathbf{V}$ correspondant aux $k$ plus grandes valeurs singulières).

Cette projection préserve au mieux la variance des données: les $k$ composantes principales capturent la direction où les données varient le plus. La reconstruction à partir de cette représentation compressée s’écrit:

L’erreur de reconstruction est minimale parmi toutes les projections linéaires sur un sous-espace de dimension $k$.

Lien entre Ridge et ACP: Les deux approches traitent le même problème (les directions à faible valeur singulière sont bruitées) mais différemment:

Ridge est appropriée pour la régression supervisée, où même les petites directions peuvent contenir du signal utile pour prédire $\mathbf{y}$. L’ACP est appropriée pour la réduction de dimension non supervisée, où nous voulons une représentation compacte des données elles-mêmes.

Pourquoi $\lambda \mathbf{I}$ aide¶

Ce terme diagonal a plusieurs effets bénéfiques:

1. Amélioration du conditionnement: La matrice $\mathbf{X}^\top \mathbf{X}$ peut être mal conditionnée (ses valeurs propres varient sur plusieurs ordres de grandeur) ou même singulière. L’ajout de $\lambda \mathbf{I}$ augmente toutes les valeurs propres de $\lambda$, rendant la matrice inversible et mieux conditionnée.

2. Rétrécissement des coefficients (shrinkage): Comme nous l’avons vu dans la section SVD ci-dessus, la solution Ridge s’écrit:

Le facteur de rétrécissement $\frac{d_j^2}{d_j^2 + \lambda}$ est toujours inférieur à 1, ce qui “rétrécit” chaque composante vers zéro. L’effet est différencié selon les directions:

• Pour une grande valeur singulière $d_j$ (fort signal), le facteur $\frac{d_j^2}{d_j^2 + \lambda}$ reste proche de 1 même pour des valeurs modérées de $\lambda$. La direction est peu affectée.

• Pour une petite valeur singulière $d_j$ (faible signal), le facteur $\frac{d_j^2}{d_j^2 + \lambda}$ décroît rapidement avec $\lambda$. La direction est fortement pénalisée.

Pour visualiser ce rétrécissement et comprendre son effet, examinons un exemple concret. L’animation suivante montre simultanément trois perspectives sur la régularisation Ridge: les données et la droite ajustée, le paysage de perte avec la contrainte, et les facteurs de rétrécissement.

Source

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.patches import Circle
from IPython.display import Image

# Générer des données de régression simple
np.random.seed(42)
n = 30

# Une seule caractéristique pour visualisation claire
x = np.random.uniform(-2, 2, n)
# Relation linéaire avec bruit
theta_true = 1.5
y = theta_true * x + np.random.normal(0, 0.8, n)

# Ajouter une caractéristique corrélée (pour créer de la colinéarité)
x2 = 0.9 * x + 0.3 * np.random.randn(n)

# Matrice de design avec les deux caractéristiques
X = np.column_stack([x, x2])

# Solution MCO
theta_ols = np.linalg.lstsq(X, y, rcond=None)[0]

# SVD pour analyse
U, d_svd, Vt = np.linalg.svd(X, full_matrices=False)
V = Vt.T

# Fonction pour calculer la solution Ridge
def ridge_solution(X, y, lam):
    n_features = X.shape[1]
    return np.linalg.solve(X.T @ X + lam * np.eye(n_features), X.T @ y)

# Préparer la grille pour les contours RSS
theta1_range = np.linspace(-0.5, 3, 100)
theta2_range = np.linspace(-1.5, 2, 100)
T1, T2 = np.meshgrid(theta1_range, theta2_range)

# Calculer RSS pour chaque point de la grille
RSS = np.zeros_like(T1)
for i in range(T1.shape[0]):
    for j in range(T1.shape[1]):
        theta = np.array([T1[i, j], T2[i, j]])
        residuals = y - X @ theta
        RSS[i, j] = np.sum(residuals**2)

# Créer la figure avec trois panneaux
fig = plt.figure(figsize=(15, 5))

# === Panneau 1: Données et droite ajustée ===
ax1 = fig.add_subplot(1, 3, 1)

# Données
ax1.scatter(x, y, c='tab:blue', s=50, alpha=0.7, label='Données', zorder=3)

# Grille pour tracer les droites
x_grid = np.linspace(-2.5, 2.5, 100)

# Droite MCO (fixe) - on utilise seulement theta1 car x et x2 sont très corrélés
# La prédiction effective est environ (theta1 + 0.9*theta2) * x
slope_ols = theta_ols[0] + 0.9 * theta_ols[1]  # Pente effective
y_ols = slope_ols * x_grid
ax1.plot(x_grid, y_ols, 'k-', linewidth=2, alpha=0.7, label='MCO')

# Droite Ridge (animée)
line_ridge, = ax1.plot([], [], '-', color='tab:orange', linewidth=2.5, label='Ridge')

# Ligne horizontale (prédiction = moyenne, lambda infini)
y_mean = np.mean(y)
ax1.axhline(y_mean, color='gray', linestyle=':', alpha=0.5, label=f'Moyenne ($\\lambda \\to \\infty$)')

ax1.set_xlabel('$x$')
ax1.set_ylabel('$y$')
ax1.set_title('Données et droite de régression')
ax1.legend(loc='upper left', fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-2.5, 2.5)
ax1.set_ylim(-4, 5)

# Texte pour les coefficients
coef_text = ax1.text(0.98, 0.02, '', transform=ax1.transAxes, fontsize=10, 
                     ha='right', va='bottom',
                     bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

# === Panneau 2: Paysage de perte ===
ax2 = fig.add_subplot(1, 3, 2)

# Contours RSS (ellipses centrées sur OLS)
levels = np.percentile(RSS.flatten(), [5, 15, 30, 50, 70, 85, 95])
contours = ax2.contour(T1, T2, RSS, levels=levels, colors='gray', alpha=0.6)
ax2.clabel(contours, inline=True, fontsize=8, fmt='%.0f')

# Solution MCO (fixe)
ax2.plot(theta_ols[0], theta_ols[1], 'ko', markersize=12, label='MCO', zorder=5)
ax2.annotate('MCO', xy=(theta_ols[0], theta_ols[1]), 
             xytext=(theta_ols[0] + 0.2, theta_ols[1] + 0.2),
             fontsize=11, ha='left')

# Origine = coefficients nuls (prédiction constante)
ax2.plot(0, 0, 'k+', markersize=15, markeredgewidth=2, zorder=4)
ax2.annotate('$\\boldsymbol{\\theta} = 0$\n(pente nulle)', xy=(0, 0), 
             xytext=(-0.4, -1.2), fontsize=9, ha='center', color='gray')

# Cercle de contrainte Ridge (animé)
circle_ridge = Circle((0, 0), radius=np.linalg.norm(theta_ols), 
                       fill=False, edgecolor='tab:orange', linewidth=2.5, 
                       linestyle='-', alpha=0.8, zorder=3)
ax2.add_patch(circle_ridge)

# Solution Ridge (animée)
point_ridge, = ax2.plot([], [], 'o', color='tab:orange', markersize=10, 
                        label='Ridge', zorder=6)

# Chemin de régularisation
lambda_path = np.logspace(-3, 1.5, 50)
theta_path = np.array([ridge_solution(X, y, l) for l in lambda_path])
ax2.plot(theta_path[:, 0], theta_path[:, 1], 'tab:orange', linewidth=1.5, 
         alpha=0.4, linestyle='--', label='Chemin')

ax2.set_xlabel('$\\theta_1$')
ax2.set_ylabel('$\\theta_2$')
ax2.set_title('Paysage RSS et chemin de régularisation')
ax2.legend(loc='upper right', fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-0.5, 3)
ax2.set_ylim(-1.5, 2)
ax2.set_aspect('equal')

# Texte pour lambda
lambda_text = ax2.text(0.02, 0.98, '', transform=ax2.transAxes, fontsize=11, 
                       va='top', ha='left',
                       bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))

# === Panneau 3: Facteurs de rétrécissement ===
ax3 = fig.add_subplot(1, 3, 3)

lambda_range = np.linspace(0, 10, 200)
shrink1_curve = d_svd[0]**2 / (d_svd[0]**2 + lambda_range)
shrink2_curve = d_svd[1]**2 / (d_svd[1]**2 + lambda_range)

ax3.plot(lambda_range, shrink1_curve, 'b-', linewidth=2, 
         label=f'Direction forte ($d_1={d_svd[0]:.1f}$)')
ax3.plot(lambda_range, shrink2_curve, 'r-', linewidth=2, 
         label=f'Direction faible ($d_2={d_svd[1]:.2f}$)')

# Zone de surapprentissage et sous-apprentissage
ax3.axvspan(0, 0.5, alpha=0.1, color='red', label='Surapprentissage')
ax3.axvspan(5, 10, alpha=0.1, color='blue', label='Sous-apprentissage')

point_shrink1, = ax3.plot([], [], 'bo', markersize=10, zorder=3)
point_shrink2, = ax3.plot([], [], 'ro', markersize=10, zorder=3)

ax3.axhline(1.0, color='gray', linestyle='--', alpha=0.5)
ax3.set_xlabel('$\\lambda$')
ax3.set_ylabel('Facteur de rétrécissement')
ax3.set_title('Rétrécissement par direction SVD')
ax3.legend(loc='center right', fontsize=8)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 10)
ax3.set_ylim(0, 1.1)

shrink_text = ax3.text(0.02, 0.5, '', transform=ax3.transAxes, fontsize=10, 
                       va='center', ha='left',
                       bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()

# Fonction d'animation
def animate(frame):
    if frame < 80:
        lam = (frame / 80) * 10
    else:
        lam = 10
    
    # Solution Ridge
    theta_ridge = ridge_solution(X, y, lam)
    
    # Panneau 1: Mettre à jour la droite
    slope_ridge = theta_ridge[0] + 0.9 * theta_ridge[1]
    y_ridge = slope_ridge * x_grid
    line_ridge.set_data(x_grid, y_ridge)
    coef_text.set_text(f'Pente MCO: {slope_ols:.2f}\nPente Ridge: {slope_ridge:.2f}')
    
    # Panneau 2: Mettre à jour le cercle et le point
    norm_ridge = np.linalg.norm(theta_ridge)
    circle_ridge.set_radius(norm_ridge)
    point_ridge.set_data([theta_ridge[0]], [theta_ridge[1]])
    lambda_text.set_text(f'$\\lambda = {lam:.1f}$')
    
    # Panneau 3: Mettre à jour les points de rétrécissement
    shrink1 = d_svd[0]**2 / (d_svd[0]**2 + lam)
    shrink2 = d_svd[1]**2 / (d_svd[1]**2 + lam)
    point_shrink1.set_data([lam], [shrink1])
    point_shrink2.set_data([lam], [shrink2])
    shrink_text.set_text(f'Facteur dir. 1: {shrink1:.2f}\nFacteur dir. 2: {shrink2:.2f}')
    
    return (line_ridge, point_ridge, circle_ridge, lambda_text, 
            point_shrink1, point_shrink2, coef_text, shrink_text)

# Créer l'animation
anim = FuncAnimation(fig, animate, frames=90, interval=80, blit=False, repeat=True)
anim.save('_static/ridge_geometry.gif', writer='pillow', fps=12, dpi=100)
plt.close()

# Afficher le GIF
Image(filename='_static/ridge_geometry.gif')

L’animation relie trois perspectives sur la régularisation Ridge lorsque $\lambda$ augmente de 0 à 10:

Panneau de gauche (données et ajustement): Les points bleus sont les données d’entraînement. La droite noire est l’ajustement MCO ($\lambda = 0$), la droite orange est l’ajustement Ridge. À mesure que $\lambda$ augmente, la pente de la droite Ridge diminue, se rapprochant de la ligne horizontale (prédiction constante égale à la moyenne). C’est le rétrécissement vers zéro: Ridge “tire” les coefficients vers l’origine, ce qui réduit la pente.

Panneau central (paysage de perte): Chaque point de ce plan représente un choix de coefficients $(\theta_1, \theta_2)$. Les contours gris montrent la fonction de coût RSS: plus on est proche du point noir (MCO), plus l’erreur sur les données d’entraînement est faible. L’ellipse est allongée car $x_1$ et $x_2$ sont corrélées (colinéarité). L’origine $\boldsymbol{\theta} = (0, 0)$ correspond à une pente nulle (prédiction constante). Le cercle orange montre la norme de la solution Ridge courante $\|\hat{\boldsymbol{\theta}}_{\text{ridge}}\|_2$. La formulation pénalisée $\text{RSS} + \lambda\|\boldsymbol{\theta}\|^2$ est équivalente à la formulation contrainte $\min \text{RSS}$ sous $\|\boldsymbol{\theta}\|^2 \leq t$, où $\lambda$ joue le rôle du multiplicateur de Lagrange: pour chaque $\lambda$, il existe un $t$ tel que les deux problèmes ont la même solution. À mesure que $\lambda$ augmente, la solution se déplace le long du chemin de régularisation vers l’origine.

Panneau de droite (rétrécissement différencié): La direction “forte” (grande valeur singulière $d_1$, où les données sont dispersées) est peu affectée par la régularisation. La direction “faible” (petite valeur singulière $d_2$, direction de colinéarité) est rétrécit beaucoup plus rapidement. C’est le cœur de l’effet Ridge: pénaliser davantage les directions où le signal est faible et l’estimation instable.

L’intuition géométrique est la suivante: quand les données sont colinéaires, l’ellipse RSS est très allongée. De petites perturbations dans les données causent de grands déplacements de la solution MCO le long de l’axe allongé. La pénalité Ridge ajoute un terme $\lambda\|\boldsymbol{\theta}\|^2$ qui « tire » la solution vers l’origine. Dans la formulation contrainte équivalente, cela correspond à chercher le minimum de RSS à l’intérieur d’une boule de rayon $\sqrt{t}$. Plus la boule est petite (plus $\lambda$ est grand), plus la solution est proche de l’origine et donc plus stable.

3. Stabilité numérique: Quand $\mathbf{X}^\top \mathbf{X}$ est presque singulière, de petites perturbations dans les données causent de grandes variations dans $\hat{\boldsymbol{\theta}}_{\text{MCO}}$. La régularisation réduit cette sensibilité.

Classes de modèles et expansion de caractéristiques¶

Les exemples de régularisation ci-dessus utilisaient des caractéristiques polynomiales: au lieu de prédire $y$ directement à partir de $x$, nous avons construit des caractéristiques $[1, x, x^2, \ldots, x^{15}]$ et appliqué un modèle linéaire dans cet espace étendu. Cette technique s’appelle l’expansion de caractéristiques et mérite d’être formalisée.

Trois familles de modèles¶

Situons les modèles linéaires dans une hiérarchie plus large. Nous distinguons trois familles de complexité croissante:

1. Modèles linéaires: $f(\mathbf{x}; \boldsymbol{\theta}) = \boldsymbol{\theta}^\top \mathbf{x} + b$. La sortie est une combinaison linéaire des entrées. Simple, interprétable, mais limité aux relations linéaires.

2. Modèles à expansion de caractéristiques: $f(\mathbf{x}; \boldsymbol{\theta}) = \boldsymbol{\theta}^\top \boldsymbol{\phi}(\mathbf{x}) + b$, où $\boldsymbol{\phi}: \mathbb{R}^d \to \mathbb{R}^D$ est une transformation non linéaire fixée à l’avance (par exemple, polynomiale). Le modèle reste linéaire dans les paramètres $\boldsymbol{\theta}$, ce qui facilite l’optimisation, mais peut capturer des relations non linéaires en $\mathbf{x}$. L’espace de redescription a souvent une dimension $D \gg d$.

3. Réseaux de neurones: $f(\mathbf{x}; \boldsymbol{\theta}) = f_K(f_{K-1}(\cdots f_1(\mathbf{x}; \boldsymbol{\theta}_1); \boldsymbol{\theta}_{K-1}); \boldsymbol{\theta}_K)$. Une composition de $K$ fonctions non linéaires, chacune avec ses propres paramètres. Contrairement aux modèles à expansion fixe, les réseaux de neurones apprennent la représentation $\boldsymbol{\phi}$ en même temps que les paramètres $\boldsymbol{\theta}$.

Cette progression capture l’évolution historique du domaine: des modèles linéaires classiques aux méthodes à noyaux (expansion implicite), puis aux réseaux profonds qui apprennent leurs propres représentations. Nous verrons les réseaux de neurones en détail dans les chapitres suivants; concentrons-nous ici sur les deux premières familles.

Expansion de caractéristiques¶

Pour capturer des relations non linéaires tout en gardant un modèle linéaire dans les paramètres, nous transformons les entrées. En régression polynomiale, nous appliquons une fonction $\phi: \mathbb{R} \to \mathbb{R}^{k+1}$:

La prédiction devient $f(x; \boldsymbol{\theta}) = \boldsymbol{\theta}^\top \phi(x)$. Le modèle est polynomial en $x$ mais linéaire en $\boldsymbol{\theta}$, ce qui permet d’utiliser les mêmes algorithmes d’optimisation (MCO, Ridge).

Le degré $k$ contrôle la capacité du modèle. Avec $k = 1$, nous avons une droite. Avec $k$ élevé, le polynôme peut osciller pour passer par tous les points d’entraînement. Avec $k = N - 1$, nous pouvons interpoler exactement les $N$ points: le risque empirique atteint zéro. Mais un polynôme qui passe exactement par les points d’entraînement n’a aucune raison de bien prédire les nouveaux points.

Illustrons ce phénomène avec les données de freinage. Nous ajustons des polynômes de degrés 1, 2, 5 et 15, et comparons leurs erreurs sur les ensembles d’entraînement et de test.

Source

import numpy as np
import matplotlib.pyplot as plt
import warnings

# Suppress polyfit warnings for high-degree polynomials (expected for this demo)
warnings.filterwarnings('ignore', message='Polyfit may be poorly conditioned')

# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
                  14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
                  20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
                 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
                 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)

# Train/test split
np.random.seed(42)
indices = np.random.permutation(len(speed))
train_idx, test_idx = indices[:35], indices[35:]
speed_train, dist_train = speed[train_idx], dist[train_idx]
speed_test, dist_test = speed[test_idx], dist[test_idx]

degrees_to_plot = [1, 2, 5, 15]
degrees_eval = range(1, 16)
fig, axes = plt.subplots(2, 2, figsize=(10, 8))

# Pre-compute all errors for the summary plot later
all_train_errors = []
all_test_errors = []
for deg in degrees_eval:
    coeffs = np.polyfit(speed_train, dist_train, deg)
    all_train_errors.append(np.mean((dist_train - np.polyval(coeffs, speed_train))**2))
    all_test_errors.append(np.mean((dist_test - np.polyval(coeffs, speed_test))**2))

for ax, deg in zip(axes.flat, degrees_to_plot):
    # Fit polynomial
    coeffs = np.polyfit(speed_train, dist_train, deg)
    
    # Predictions
    pred_train = np.polyval(coeffs, speed_train)
    pred_test = np.polyval(coeffs, speed_test)
    
    # MSE
    mse_train = np.mean((dist_train - pred_train)**2)
    mse_test = np.mean((dist_test - pred_test)**2)
    
    # Plot
    ax.scatter(speed_train, dist_train, alpha=0.6, s=30, label='Entraînement')
    ax.scatter(speed_test, dist_test, alpha=0.6, s=30, marker='s', label='Test')
    
    speed_grid = np.linspace(3, 26, 200)
    pred_grid = np.polyval(coeffs, speed_grid)
    # Clip extreme predictions for visualization
    pred_grid = np.clip(pred_grid, -50, 200)
    ax.plot(speed_grid, pred_grid, 'k-', alpha=0.7)
    
    ax.set_xlim(3, 26)
    ax.set_ylim(-20, 150)
    ax.set_xlabel('Vitesse (mph)')
    ax.set_ylabel('Distance (ft)')
    ax.set_title(f'Degré {deg}: Entr. EQM={mse_train:.1f}, Test EQM={mse_test:.1f}')
    if deg == 1:
        ax.legend()

plt.tight_layout()

Le polynôme de degré 1 (droite) ne capture pas la courbure des données: c’est du sous-apprentissage. Le polynôme de degré 2 capture bien la relation quadratique. Le polynôme de degré 5 commence à osciller. Le polynôme de degré 15 passe près de tous les points d’entraînement, mais ses oscillations produisent des prédictions absurdes entre les points: c’est du surapprentissage.

Source

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(degrees_eval, all_train_errors, 'o-', linewidth=2, label='Erreur entraînement')
ax.plot(degrees_eval, all_test_errors, 's-', linewidth=2, label='Erreur test')

# Utiliser une échelle logarithmique car l'erreur de test explose
ax.set_yscale('log')

ax.set_xlabel('Degré du polynôme (complexité)')
ax.set_ylabel('EQM (échelle log)')
ax.set_xticks(range(1, 16, 2))
ax.grid(True, which="both", ls="-", alpha=0.2)
ax.legend()

ax.set_title('Compromis biais-variance')
plt.tight_layout()

L’erreur d’entraînement diminue avec le degré du polynôme. L’erreur de test diminue d’abord (quand le modèle gagne en expressivité), puis augmente (quand le modèle commence à mémoriser le bruit). Le meilleur modèle se trouve à l’intersection de ces deux tendances. La régularisation Ridge, vue précédemment, est une alternative au choix du degré: elle permet d’utiliser un modèle de haute capacité tout en contrôlant le surapprentissage.

Décomposition biais-variance¶

Ce compromis peut être formalisé mathématiquement. Supposons que les données suivent le modèle $y = f^*(\mathbf{x}) + \epsilon$, où $f^*$ est la vraie fonction (le prédicteur de Bayes optimal pour la perte quadratique), et $\epsilon$ est un bruit de moyenne nulle et de variance $\sigma^2$.

Clarification sur les variables aléatoires: Dans cette analyse, il y a trois sources d’aléa qu’il faut bien distinguer:

• $f^*(\mathbf{x})$ est une valeur déterministe (fixe) pour chaque $\mathbf{x}$. C’est la vraie fonction sous-jacente, qui ne dépend d’aucun tirage aléatoire.

• $\epsilon$ est une variable aléatoire représentant le bruit d’observation. C’est l’aléa intrinsèque aux données.

• $\hat{f}$ est notre estimateur, une fonction apprise à partir de l’échantillon d’entraînement $\mathcal{D}$. Comme $\mathcal{D}$ est tiré aléatoirement, différents tirages de $\mathcal{D}$ produisent différentes fonctions $\hat{f}$. Pour un point test $\mathbf{x}$ fixé, la prédiction $\hat{f}(\mathbf{x})$ est donc une variable aléatoire (un scalaire qui varie selon $\mathcal{D}$), même si $\hat{f}$ elle-même est une fonction.

Autrement dit, l’espérance $\mathbb{E}[\hat{f}(\mathbf{x})]$ moyenne sur tous les échantillons d’entraînement possibles: si nous pouvions répéter l’expérience d’apprentissage un grand nombre de fois avec différents $\mathcal{D}$, quelle serait la prédiction moyenne pour ce $\mathbf{x}$?

L’erreur quadratique moyenne, moyennée sur les échantillons d’entraînement possibles et le bruit, se décompose comme suit:

En développant le carré et utilisant $\mathbb{E}[\epsilon] = 0$ ainsi que l’indépendance entre $\epsilon$ et $\hat{f}$:

Pour le premier terme, ajoutons et retranchons $\mathbb{E}[\hat{f}(\mathbf{x})]$:

En développant et utilisant $\mathbb{E}[\hat{f}(\mathbf{x}) - \mathbb{E}[\hat{f}(\mathbf{x})]] = 0$:

Nous obtenons la décomposition biais-variance:

Chaque terme a une interprétation précise:

• Biais²: L’écart entre la prédiction moyenne de notre estimateur et la vraie fonction $f^*$. Un modèle trop simple (classe $\mathcal{H}$ trop restrictive) aura un biais élevé car il ne peut pas approcher $f^*$.

• Variance: La sensibilité de notre estimateur à l’échantillon d’entraînement particulier. Un modèle trop complexe aura une variance élevée car de petites variations dans les données causent de grandes variations dans les prédictions.

• $\sigma^2$: Le bruit irréductible, inhérent aux données. Aucun estimateur ne peut faire mieux que cette erreur.

Cette décomposition explique pourquoi l’erreur de test a une forme en U: à faible complexité, le biais domine; à haute complexité, la variance domine. Le minimum se trouve au point où la somme des deux est minimale.

Intuition géométrique: pourquoi la dimension supérieure aide¶

L’expansion de caractéristiques semble être un simple changement de variables, mais elle cache une idée géométrique profonde. Pour comprendre pourquoi projeter les données dans un espace de dimension supérieure permet de capturer des relations non linéaires, examinons d’abord le cas de la régression, puis celui de la classification.

Le plan caché derrière la parabole¶

Considérons une régression quadratique: $f(x) = \theta_0 + \theta_1 x + \theta_2 x^2$. Cette fonction est non linéaire en $x$ (c’est une parabole), mais linéaire dans les paramètres $\boldsymbol{\theta} = (\theta_0, \theta_1, \theta_2)$. Que signifie cette distinction géométriquement?

Introduisons l’espace des caractéristiques $\phi(x) = (1, x, x^2)$. Chaque valeur de $x$ correspond à un point dans $\mathbb{R}^3$. Ces points ne sont pas dispersés arbitrairement: ils vivent sur une courbe particulière, la courbe des moments (moment curve), qui ressemble à une rampe tordue.

Source

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Créer la courbe des moments: (x, x², x³) pour visualisation
# On utilise (1, x, x²) mais on projette sur (x, x², y) pour la visualisation
fig = plt.figure(figsize=(12, 5))

# Gauche: les fonctions de base
ax1 = fig.add_subplot(121)
x = np.linspace(-2, 2, 100)
ax1.plot(x, np.ones_like(x), 'b-', linewidth=2, label=r'$\phi_0(x) = 1$')
ax1.plot(x, x, 'orange', linewidth=2, label=r'$\phi_1(x) = x$')
ax1.plot(x, x**2, 'g-', linewidth=2, label=r'$\phi_2(x) = x^2$')
ax1.axhline(0, color='gray', linewidth=0.5)
ax1.axvline(0, color='gray', linewidth=0.5)
ax1.set_xlabel('$x$')
ax1.set_ylabel(r'$\phi_j(x)$')
ax1.set_title('Les fonctions de base')
ax1.legend()
ax1.set_ylim(-2.5, 4.5)
ax1.grid(True, alpha=0.3)

# Droite: combinaisons linéaires
ax2 = fig.add_subplot(122)
x = np.linspace(-2, 2, 100)

# Différentes combinaisons de coefficients
combinations = [
    ((1, 0, 0), 'Constante: $1$'),
    ((0, 1, 0), 'Linéaire: $x$'),
    ((0, 0, 1), 'Quadratique: $x^2$'),
    ((1, -0.5, 0.5), 'Combinaison: $1 - 0.5x + 0.5x^2$'),
]

colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red']
for (theta0, theta1, theta2), label in combinations:
    y = theta0 + theta1 * x + theta2 * x**2
    ax2.plot(x, y, linewidth=2, label=label)

ax2.axhline(0, color='gray', linewidth=0.5)
ax2.axvline(0, color='gray', linewidth=0.5)
ax2.set_xlabel('$x$')
ax2.set_ylabel('$f(x)$')
ax2.set_title('Combinaisons linéaires des fonctions de base')
ax2.legend(loc='upper center')
ax2.set_ylim(-2.5, 4.5)
ax2.grid(True, alpha=0.3)

plt.tight_layout()

Les fonctions de base $\{1, x, x^2\}$ sont les “ingrédients” du modèle. La régression polynomiale cherche les coefficients $\theta_0, \theta_1, \theta_2$ qui mélangent ces ingrédients de façon optimale. Chaque combinaison produit une courbe différente, mais toutes sont des paraboles (ou des cas dégénérés: droites, constantes).

Voici l’insight géométrique clé: dans l’espace $(x, x^2, y)$, le modèle $y = \theta_0 + \theta_1 x + \theta_2 x^2$ définit un plan. La parabole que nous voyons dans le graphique $(x, y)$ est simplement la projection de ce plan sur notre espace de visualisation.

Source

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from IPython.display import Image

# Générer des données avec relation quadratique
np.random.seed(42)
n_points = 40
x_data = np.random.uniform(-1.8, 1.8, n_points)
y_true = 0.5 + 0.3 * x_data + 0.8 * x_data**2
y_data = y_true + np.random.normal(0, 0.3, n_points)

# Ajuster le modèle quadratique
X_design = np.column_stack([np.ones(n_points), x_data, x_data**2])
theta = np.linalg.lstsq(X_design, y_data, rcond=None)[0]

# Créer l'animation
fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(111, projection='3d')

# Grille pour le plan de régression
x_grid = np.linspace(-2, 2, 20)
x2_grid = np.linspace(0, 4, 20)
X_plane, X2_plane = np.meshgrid(x_grid, x2_grid)
Y_plane = theta[0] + theta[1] * X_plane + theta[2] * X2_plane

# Surface parabolique z = x²
x_surf = np.linspace(-2, 2, 30)
X_surf, Y_surf_temp = np.meshgrid(x_surf, np.linspace(-1, 5, 30))
Z_surf = X_surf**2

def init():
    ax.clear()
    return []

def animate(frame):
    ax.clear()
    
    # Animation en trois phases
    if frame < 15:
        # Phase 1: Vue 2D (de côté, cachant x²)
        elev = 0
        azim = 0
        show_plane = False
        show_surface = False
        title = 'Vue 2D: régression quadratique'
    elif frame < 35:
        # Phase 2: Rotation révélant la 3ème dimension
        progress = (frame - 15) / 20
        elev = progress * 25
        azim = progress * 45
        show_plane = False
        show_surface = True
        title = 'Rotation: découverte de la dimension $x^2$...'
    else:
        # Phase 3: Vue 3D avec le plan
        elev = 25
        azim = 45 + (frame - 35) * 2
        show_plane = True
        show_surface = True
        title = r'Le plan $y = \theta_0 + \theta_1 x + \theta_2 x^2$ dans lespace 3D'
    
    ax.view_init(elev=elev, azim=azim)
    
    # Surface parabolique (rampe z = x²)
    if show_surface:
        ax.plot_surface(X_surf, Z_surf, Y_surf_temp, alpha=0.1, color='gray')
    
    # Points de données dans l'espace 3D: (x, x², y)
    ax.scatter(x_data, x_data**2, y_data, c='tab:blue', s=50, alpha=0.8, 
               label='Données', depthshade=True)
    
    # Plan de régression
    if show_plane:
        ax.plot_surface(X_plane, X2_plane, Y_plane, alpha=0.4, color='tab:orange',
                       label='Plan de régression')
    
    # Courbe de régression sur la surface z = x²
    x_curve = np.linspace(-1.8, 1.8, 100)
    y_curve = theta[0] + theta[1] * x_curve + theta[2] * x_curve**2
    ax.plot(x_curve, x_curve**2, y_curve, 'r-', linewidth=3, 
            label='Courbe ajustée')
    
    ax.set_xlabel('$x$')
    ax.set_ylabel('$x^2$')
    ax.set_zlabel('$y$')
    ax.set_xlim(-2, 2)
    ax.set_ylim(0, 4)
    ax.set_zlim(-1, 5)
    ax.set_title(title, fontsize=11)
    
    return []

anim = FuncAnimation(fig, animate, init_func=init, frames=70, interval=100, blit=True)
anim.save('_static/regression_plane_3d.gif', writer='pillow', fps=10, dpi=100)
plt.close()

# Afficher le GIF
Image(filename='_static/regression_plane_3d.gif')

L’animation révèle la structure cachée de la régression polynomiale:

1. Vue 2D initiale: On voit les données et la parabole ajustée, comme dans un graphique classique.

2. Rotation: En faisant pivoter la vue, on découvre que chaque point $(x_i, y_i)$ vit en réalité dans un espace 3D aux coordonnées $(x_i, x_i^2, y_i)$. La surface grise représente la “rampe” $z = x^2$ sur laquelle tous les points sont contraints de vivre.

3. Le plan de régression: Le modèle $y = \theta_0 + \theta_1 x + \theta_2 x^2$ est un plan (en orange) dans cet espace 3D. Ce plan est choisi pour minimiser les distances verticales aux points.

4. La courbe ajustée: La parabole rouge est l’intersection du plan avec la rampe $z = x^2$. C’est ce que nous voyons quand nous projetons le tout sur le plan $(x, y)$.

Cette perspective unifie le “linéaire dans les paramètres” et le “non linéaire en $x$”:

• Linéaire dans les paramètres: Le modèle est un plan (objet linéaire) dans l’espace des caractéristiques.

• Non linéaire en $x$: La contrainte $z = x^2$ force les données à vivre sur une surface courbe, et l’intersection du plan avec cette surface produit une courbe.

De la régression à la classification¶

La même intuition s’applique en classification, avec une différence: au lieu de chercher un plan qui ajuste les données, on cherche un hyperplan qui les sépare. Voyons comment l’expansion de caractéristiques transforme des données non linéairement séparables en données linéairement séparables.

De 1D à 2D: séparer l’inséparable¶

Considérons des points sur une droite, répartis en deux classes: les points bleus au centre, les points orange aux extrémités. Aucun seuil unique ne peut séparer ces deux classes.

Source

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# Classe bleue: points au centre
x_blue = np.random.uniform(-0.5, 0.5, 15)
# Classe orange: points aux extrémités
x_orange = np.concatenate([np.random.uniform(-1.5, -0.8, 8), 
                           np.random.uniform(0.8, 1.5, 8)])

fig, ax = plt.subplots(figsize=(10, 2))
ax.scatter(x_blue, np.zeros_like(x_blue), c='tab:blue', s=80, label='Classe A', zorder=3)
ax.scatter(x_orange, np.zeros_like(x_orange), c='tab:orange', s=80, label='Classe B', zorder=3)
ax.axhline(0, color='gray', linewidth=0.5, zorder=1)
ax.set_xlim(-2, 2)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('$x$')
ax.set_yticks([])
ax.legend(loc='upper right')
ax.set_title('Données 1D: aucun seuil ne sépare les deux classes')
plt.tight_layout()

Appliquons maintenant l’expansion $\phi(x) = (x, x^2)$. Chaque point est projeté sur une parabole dans l’espace 2D. Les points proches de zéro (classe bleue) restent bas, tandis que les points éloignés de zéro (classe orange) montent.

Source

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Gauche: espace original avec tentative de séparation
ax = axes[0]
ax.scatter(x_blue, np.zeros_like(x_blue), c='tab:blue', s=80, zorder=3)
ax.scatter(x_orange, np.zeros_like(x_orange), c='tab:orange', s=80, zorder=3)
ax.axhline(0, color='gray', linewidth=0.5, zorder=1)
ax.axvline(0.3, color='red', linestyle='--', linewidth=2, label='Seuil?')
ax.set_xlim(-2, 2)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('$x$')
ax.set_yticks([])
ax.set_title('Espace original: pas de séparation linéaire')
ax.legend()

# Droite: espace transformé
ax = axes[1]
ax.scatter(x_blue, x_blue**2, c='tab:blue', s=80, zorder=3, label='Classe A')
ax.scatter(x_orange, x_orange**2, c='tab:orange', s=80, zorder=3, label='Classe B')

# Ligne de séparation dans l'espace transformé
x_line = np.linspace(-2, 2, 100)
threshold = 0.6
ax.axhline(threshold, color='green', linestyle='-', linewidth=2, label='Frontière linéaire')
ax.fill_between(x_line, 0, threshold, alpha=0.1, color='blue')
ax.fill_between(x_line, threshold, 2.5, alpha=0.1, color='orange')

# Parabole de référence
x_curve = np.linspace(-1.6, 1.6, 100)
ax.plot(x_curve, x_curve**2, 'k--', alpha=0.3, linewidth=1)

ax.set_xlim(-2, 2)
ax.set_ylim(-0.1, 2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$x^2$')
ax.set_title(r'Espace transformé $\phi(x) = (x, x^2)$: séparation linéaire!')
ax.legend()

plt.tight_layout()

Dans l’espace transformé, une simple droite horizontale sépare les deux classes. Cette droite correspond, dans l’espace original, à deux seuils: $x^2 < 0.6$, soit $|x| < \sqrt{0.6} \approx 0.77$. L’expansion de caractéristiques a transformé une frontière de décision non linéaire (un intervalle) en une frontière linéaire (une droite).

De 2D à 3D: soulever pour séparer¶

Passons à un exemple plus visuel. Considérons deux classes disposées en cercles concentriques: la classe bleue forme un disque central, la classe orange forme un anneau extérieur. Aucune droite ne peut séparer ces deux régions.

Source

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

# Classe bleue: disque central
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)

# Classe orange: anneau extérieur
n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(x_blue, y_blue, c='tab:blue', s=40, alpha=0.7, label='Classe A')
ax.scatter(x_orange, y_orange, c='tab:orange', s=40, alpha=0.7, label='Classe B')

# Montrer qu'une droite ne peut pas séparer
theta_line = np.pi/4
x_line = np.linspace(-2, 2, 100)
y_line = np.tan(theta_line) * x_line
ax.plot(x_line, y_line, 'r--', linewidth=2, alpha=0.7, label='Droite?')

ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_aspect('equal')
ax.legend()
ax.set_title('Cercles concentriques: pas de séparation linéaire en 2D')
plt.tight_layout()

Appliquons l’expansion $\phi(x_1, x_2) = (x_1, x_2, x_1^2 + x_2^2)$. La troisième coordonnée est le carré de la distance à l’origine: $z = r^2$. Les points proches du centre (petit $r$) sont “soulevés” moins haut que les points éloignés (grand $r$).

L’animation suivante montre cette transformation. En faisant pivoter la vue, on voit que les données, une fois projetées dans l’espace 3D, deviennent séparables par un plan horizontal.

Source

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import Image

np.random.seed(123)

# Classe bleue: disque central
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)
z_blue = x_blue**2 + y_blue**2

# Classe orange: anneau extérieur  
n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)
z_orange = x_orange**2 + y_orange**2

# Créer la figure
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

# Plan de séparation
z_sep = 0.75
xx, yy = np.meshgrid(np.linspace(-1.8, 1.8, 10), np.linspace(-1.8, 1.8, 10))
zz = np.ones_like(xx) * z_sep

def init():
    ax.clear()
    return []

def animate(frame):
    ax.clear()
    
    # Élévation: commence à 90° (vue de dessus), descend à 25°
    if frame < 20:
        elev = 90 - frame * 3.25  # 90 -> 25
        azim = 45
    else:
        elev = 25
        azim = 45 + (frame - 20) * 4  # rotation horizontale
    
    ax.view_init(elev=elev, azim=azim)
    
    # Points
    ax.scatter(x_blue, y_blue, z_blue, c='tab:blue', s=30, alpha=0.8, label='Classe A')
    ax.scatter(x_orange, y_orange, z_orange, c='tab:orange', s=30, alpha=0.8, label='Classe B')
    
    # Plan de séparation (apparaît après la descente)
    if frame >= 15:
        alpha_plane = min(0.3, (frame - 15) * 0.06)
        ax.plot_surface(xx, yy, zz, alpha=alpha_plane, color='green')
    
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$x_1^2 + x_2^2$')
    ax.set_xlim(-1.8, 1.8)
    ax.set_ylim(-1.8, 1.8)
    ax.set_zlim(0, 2.5)
    
    # Titre dynamique
    if frame < 10:
        ax.set_title('Vue de dessus: cercles concentriques', fontsize=11)
    elif frame < 20:
        ax.set_title('Rotation: on découvre la structure 3D...', fontsize=11)
    else:
        ax.set_title(r'Espace $\phi(x_1,x_2) = (x_1, x_2, x_1^2+x_2^2)$: un plan sépare!', fontsize=11)
    
    return []

anim = FuncAnimation(fig, animate, init_func=init, frames=60, interval=100, blit=True)
anim.save('_static/feature_expansion_3d.gif', writer='pillow', fps=10, dpi=100)
plt.close()

# Afficher le GIF
Image(filename='_static/feature_expansion_3d.gif')

L’animation montre comment l’expansion de caractéristiques “déplie” la géométrie des données. Vue de dessus, la structure est celle des cercles concentriques originaux. Mais la troisième dimension $z = x_1^2 + x_2^2$ sépare verticalement les deux classes: le disque central reste bas, l’anneau extérieur monte. Un plan horizontal (en vert) suffit alors à séparer les classes.

Ce plan horizontal $z = 0.75$ correspond, dans l’espace original 2D, à un cercle de rayon $\sqrt{0.75} \approx 0.87$. La frontière de décision linéaire en 3D devient une frontière circulaire en 2D.

Source

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(123)

# Recréer les données (même seed que l'animation)
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)

n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(x_blue, y_blue, c='tab:blue', s=40, alpha=0.7, label='Classe A')
ax.scatter(x_orange, y_orange, c='tab:orange', s=40, alpha=0.7, label='Classe B')

# Cercle de décision (projection du plan z = 0.75)
theta_circle = np.linspace(0, 2*np.pi, 100)
r_decision = np.sqrt(0.75)
ax.plot(r_decision * np.cos(theta_circle), r_decision * np.sin(theta_circle), 
        'g-', linewidth=2.5, label=f'Frontière: $r = {r_decision:.2f}$')

ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_aspect('equal')
ax.legend()
ax.set_title('Frontière de décision projetée en 2D')
plt.tight_layout()

Le principe unificateur¶

Ces exemples de régression et de classification illustrent le même principe géométrique:

L’expansion de caractéristiques transforme un problème non linéaire en un problème linéaire dans un espace de dimension supérieure. Les modèles linéaires, simples à optimiser et à analyser, deviennent alors suffisants pour capturer des structures complexes.

En augmentant la dimension de l’espace de représentation, nous augmentons la capacité du modèle. Pour la classification, un résultat classique de géométrie affirme que $N$ points en position générale sont presque sûrement séparables par un hyperplan si la dimension de l’espace est au moins $N$. Pour la régression, un polynôme de degré $N-1$ peut interpoler exactement $N$ points.

Mais cette flexibilité a un coût: plus l’espace est grand, plus le modèle risque de mémoriser les particularités des données d’entraînement plutôt que d’apprendre la structure sous-jacente. C’est le compromis biais-variance, et c’est pourquoi la régularisation (vue précédemment) est si importante pour les modèles à haute capacité.

Évaluation et choix de modèle¶

En pratique, nous estimons le risque par le risque empirique sur un ensemble de test $\mathcal{D}_{\text{test}}$ disjoint de l’ensemble d’entraînement. Un troisième ensemble, l’ensemble de validation, sert à choisir parmi plusieurs modèles ou à régler des hyperparamètres. L’ensemble de test doit rester intact jusqu’à l’évaluation finale, pour fournir une estimation non biaisée.

Cette séparation est importante. Si nous utilisons l’ensemble de test pour faire des choix (quel modèle garder, quelle valeur d’hyperparamètre utiliser), l’estimation de performance sur ce même ensemble devient optimiste.

Hyperparamètres et validation¶

De nombreux modèles ont des hyperparamètres: des choix qui doivent être faits avant l’entraînement et qui ne sont pas appris à partir des données. Le degré $k$ d’un polynôme, le nombre de voisins dans les $k$ plus proches voisins, ou le coefficient de régularisation $\lambda$ sont des exemples d’hyperparamètres.

Un hyperparamètre mal choisi peut mener au surapprentissage (modèle trop complexe) ou au sous-apprentissage (modèle trop simple). La méthode standard pour choisir un hyperparamètre est la validation: on réserve une partie des données (typiquement 20%) comme ensemble de validation, on entraîne le modèle pour plusieurs valeurs de l’hyperparamètre, et on retient celle qui minimise l’erreur sur l’ensemble de validation.

Plus formellement, pour un hyperparamètre $h$, définissons le risque de validation:

où $\hat{f}_h(\mathcal{D}_{\text{train}})$ est le modèle entraîné avec l’hyperparamètre $h$. La recherche par grille consiste à évaluer ce risque pour un ensemble de valeurs candidates et à retenir:

Une fois $h^*$ choisi, on peut ré-entraîner le modèle sur l’ensemble des données (entraînement + validation) pour obtenir le modèle final.

Validation croisée¶

Quand les données sont peu nombreuses, réserver 20% pour la validation peut être coûteux. La validation croisée offre une alternative.

L’idée est de partitionner les données en $K$ blocs. Pour chaque bloc $k$, on entraîne le modèle sur les $K-1$ autres blocs et on évalue sur le bloc $k$. Le risque de validation croisée est la moyenne des $K$ évaluations:

où $\mathcal{D}_{-k}$ désigne toutes les données sauf le bloc $k$.

Source

import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 3))

K = 5
for fold in range(K):
    for k in range(K):
        if k == fold:
            ax.barh(fold, 1, left=k, color='C1', edgecolor='black', linewidth=1, label='Validation' if fold == 0 else '')
        else:
            ax.barh(fold, 1, left=k, color='C0', edgecolor='black', linewidth=1, label='Entraînement' if fold == 0 and k == 0 else '')

ax.set_yticks(range(K))
ax.set_yticklabels([f'Itération {k+1}' for k in range(K)])
ax.set_xticks(np.arange(K) + 0.5)
ax.set_xticklabels([f'Bloc {k+1}' for k in range(K)])
ax.set_xlabel('Blocs de données')
ax.legend(loc='upper right')
ax.set_title(f'Validation croisée à {K} blocs')
ax.set_xlim(0, K)

plt.tight_layout()

Le cas particulier $K = N$ (un bloc par exemple) est appelé validation croisée leave-one-out. Elle utilise au maximum les données disponibles, mais est coûteuse en calcul. En pratique, $K = 5$ ou $K = 10$ offrent un bon compromis.

Biais inductifs¶

Il n’existe pas de modèle universel qui fonctionne optimalement pour tous les problèmes. Ce résultat, connu sous le nom de théorème du no free lunch, affirme qu’un algorithme d’apprentissage qui performe bien sur une classe de problèmes performe nécessairement moins bien sur d’autres.

Tout modèle encode des biais inductifs: des hypothèses implicites ou explicites sur la structure du problème. La régression linéaire suppose que la relation entre entrées et sorties est linéaire. Les k plus proches voisins supposent que les points proches dans l’espace des entrées ont des sorties similaires. Les modèles plus complexes, comme les réseaux de neurones, encodent d’autres hypothèses sur la structure des données.

Ces hypothèses sont nécessaires pour que l’apprentissage soit possible. Sans elles, nous n’aurions aucune raison de croire que la performance sur l’échantillon d’entraînement prédit la performance sur de nouvelles données. Le choix du modèle et de ses hypothèses est une décision que l’algorithme ne peut pas prendre seul; elle requiert une connaissance du domaine.

Fonctions de perte de substitution¶

La perte 0-1 pose un problème pratique. Les méthodes d’optimisation itératives, comme la descente de gradient, requièrent que la fonction objectif soit différentiable. Or la perte 0-1 est constante par morceaux: sa dérivée est nulle presque partout et indéfinie aux points de discontinuité.

Nous contournons ce problème en utilisant des fonctions de perte de substitution: des approximations convexes et différentiables de la perte originale.

Pour la classification binaire, plutôt que de prédire directement une classe, les modèles produisent souvent un score $s = f(\mathbf{x})$ (un nombre réel). La prédiction de classe se fait ensuite en prenant le signe de ce score: si $s > 0$, on prédit la classe +1; si $s < 0$, on prédit la classe -1. La valeur absolue de $s$ mesure la confiance: plus $|s|$ est grand, plus le modèle est confiant dans sa prédiction.

Pour la classification binaire avec $y \in \{-1, +1\}$, la perte logistique est:

où $s = f(\mathbf{x})$ est le score produit par le modèle. Cette fonction est convexe et différentiable partout. Lorsque $y$ et $s$ ont le même signe (prédiction correcte avec confiance), la perte est faible. Lorsqu’ils ont des signes opposés (erreur), la perte croît linéairement avec l’amplitude de l’erreur.

La perte à charnière (hinge loss) est utilisée dans les machines à vecteurs de support:

Cette fonction est convexe mais non différentiable au point $y \cdot s = 1$. Elle est nulle lorsque la prédiction est correcte avec une marge suffisante ($y \cdot s \geq 1$), et croît linéairement sinon.

Ces deux fonctions majorent la perte 0-1: pour tout $y$ et $s$, nous avons $\ell_{0-1} \leq \ell_{\text{log}}$ et $\ell_{0-1} \leq \ell_{\text{hinge}}$. Minimiser ces substituts garantit donc un certain contrôle sur la perte originale.

Source

import numpy as np
import matplotlib.pyplot as plt

# Margin: y * s (positive = correct prediction, negative = error)
margin = np.linspace(-3, 3, 500)

# 0-1 loss: 1 if margin < 0, else 0
loss_01 = (margin < 0).astype(float)

# Logistic loss: log(1 + exp(-margin))
loss_log = np.log(1 + np.exp(-margin))

# Hinge loss: max(0, 1 - margin)
loss_hinge = np.maximum(0, 1 - margin)

fig, ax = plt.subplots(figsize=(8, 5))

ax.plot(margin, loss_01, 'k-', linewidth=2, label='Perte 0-1')
ax.plot(margin, loss_log, 'C0-', linewidth=2, label='Perte logistique')
ax.plot(margin, loss_hinge, 'C1-', linewidth=2, label='Perte à charnière')

ax.axvline(0, color='gray', linestyle=':', alpha=0.5)
ax.axhline(1, color='gray', linestyle=':', alpha=0.3)

ax.set_xlabel(r'Marge $y \cdot s$')
ax.set_ylabel('Perte')
ax.set_xlim(-3, 3)
ax.set_ylim(-0.1, 4)
ax.legend()
ax.set_title('Fonctions de perte de substitution comme bornes supérieures convexes')

# Annotate regions
ax.text(-1.5, 3.5, 'Erreur\n(prédiction incorrecte)', ha='center', fontsize=9, color='gray')
ax.text(1.5, 0.3, 'Correct\n(prédiction juste)', ha='center', fontsize=9, color='gray')

plt.tight_layout()

La figure montre les trois fonctions de perte en fonction de la marge $y \cdot s$. Une marge positive indique une prédiction correcte (le signe de $s$ correspond à $y$), une marge négative indique une erreur. La perte 0-1 est discontinue au point $y \cdot s = 0$. Les pertes logistique et à charnière sont continues et convexes, ce qui permet d’utiliser des méthodes d’optimisation par gradient. Elles majorent partout la perte 0-1.

Le cadre probabiliste¶

Jusqu’ici, nous avons choisi des fonctions de perte de manière ad hoc: la perte quadratique semble raisonnable pour la régression, la perte logistique pour la classification. Mais d’où viennent ces choix? Existe-t-il un principe unificateur?

Le cadre probabiliste offre une réponse: plutôt que de choisir une perte arbitraire, nous modélisons explicitement comment les données ont été générées. Cette section présente d’abord le cadre général de l’inférence bayésienne, puis développe deux approches concrètes: le maximum de vraisemblance (EMV) et le maximum a posteriori (MAP).

Le cadre bayésien¶

La statistique bayésienne propose un cadre général pour l’estimation de paramètres. Au lieu d’estimer un point unique, elle caractérise notre incertitude sur les paramètres par une distribution de probabilité.

Le théorème de Bayes nous dit comment mettre à jour nos croyances sur les paramètres $\boldsymbol{\theta}$ après avoir observé des données $\mathcal{D}$:

Chaque terme a un nom et un rôle précis:

• $p(\boldsymbol{\theta} | \mathcal{D})$ est la distribution a posteriori: notre croyance sur $\boldsymbol{\theta}$ après avoir vu les données

• $p(\boldsymbol{\theta})$ est la distribution a priori: notre croyance sur $\boldsymbol{\theta}$ avant d’observer les données

• $p(\mathcal{D} | \boldsymbol{\theta})$ est la vraisemblance: la probabilité des données pour un choix de paramètres

• $p(\mathcal{D}) = \int p(\boldsymbol{\theta}') p(\mathcal{D} | \boldsymbol{\theta}') d\boldsymbol{\theta}'$ est la vraisemblance marginale: une constante de normalisation

L’a priori encode notre connaissance préalable. Pour une pièce de monnaie, nous pourrions croire que $\theta$ est probablement proche de 0.5. L’a posteriori combine cette croyance avec l’évidence des données.

Le prédicteur de Bayes optimal¶

Commençons par un idéal théorique. Si nous connaissions la vraie distribution conjointe $p(\mathbf{x}, y)$, quelle fonction $f$ minimiserait le risque?

Puisque $p(\mathbf{x})$ est toujours positif, minimiser cette intégrale revient à minimiser, pour chaque $\mathbf{x}$, l’espérance conditionnelle de la perte. Le prédicteur de Bayes optimal est donc:

La réponse dépend de la fonction de perte choisie.

Pour la perte quadratique $\ell(y, \hat{y}) = (y - \hat{y})^2$, développons l’espérance:

C’est une fonction quadratique en $\hat{y}$. En dérivant et en posant la dérivée égale à zéro:

Le prédicteur optimal est la moyenne conditionnelle.

Pour la perte 0-1 en classification $\ell(y, \hat{y}) = \mathbf{1}[y \neq \hat{y}]$, l’espérance est:

Minimiser cette quantité revient à maximiser $P(y = \hat{y} | \mathbf{x})$, donc à choisir la classe la plus probable:

Le prédicteur optimal est le mode conditionnel. Chaque perte définit son propre prédicteur optimal.

Ce prédicteur est un repère théorique: aucun algorithme ne peut faire mieux, car il suppose l’accès à la vraie distribution. La différence entre le risque d’un prédicteur appris et ce risque de Bayes mesure ce que nous perdons en ne connaissant pas la vraie distribution.

Prédiction bayésienne et distribution prédictive a posteriori¶

En pratique, nous ne connaissons pas $p(y|\mathbf{x})$. Nous avons un modèle paramétrique $p(y|\mathbf{x}, \boldsymbol{\theta})$ et une distribution a posteriori $p(\boldsymbol{\theta}|\mathcal{D})$ sur les paramètres. L’approche bayésienne complète consiste à moyenner les prédictions sur tous les paramètres possibles, pondérés par leur probabilité a posteriori:

Cette distribution prédictive a posteriori intègre l’incertitude sur les paramètres. Elle ne s’engage pas sur une valeur unique de $\boldsymbol{\theta}$, mais considère toutes les valeurs plausibles.

Le problème: cette intégrale est rarement calculable analytiquement. Elle nécessite de sommer sur un espace de paramètres de grande dimension, ce qui est coûteux ou impossible en pratique. C’est pourquoi nous recourons souvent à des estimateurs ponctuels: plutôt que d’intégrer sur tous les $\boldsymbol{\theta}$, nous en choisissons un seul, comme l’EMV ou le MAP.

Utilité du modèle probabiliste¶

Si nous finissons souvent par utiliser un estimateur ponctuel, pourquoi adopter le cadre probabiliste? Plusieurs raisons:

1. Justifier la fonction de perte: La perte quadratique découle naturellement de l’hypothèse de bruit gaussien. La perte logarithmique vient du principe de maximum de vraisemblance. Le cadre probabiliste explique pourquoi ces choix sont raisonnables.

2. Quantifier l’incertitude: Au-delà de la prédiction ponctuelle $\hat{y} = f(\mathbf{x}; \hat{\boldsymbol{\theta}})$, nous pouvons donner un intervalle de prédiction. Sous un modèle gaussien, $y$ a environ 95% de chances de tomber dans $[f(\mathbf{x}) - 2\sigma, f(\mathbf{x}) + 2\sigma]$.

3. Comparer des modèles: La vraisemblance marginale $p(\mathcal{D})$ permet de comparer des modèles de complexités différentes, pénalisant automatiquement les modèles trop complexes.

4. Ouvrir la porte à l’inférence complète: Quand les ressources le permettent (méthodes de Monte Carlo, inférence variationnelle), nous pouvons approximer la distribution prédictive complète plutôt que de nous limiter à un point.

Maximum de vraisemblance¶

Le maximum de vraisemblance est la première approche concrète dans le cadre probabiliste: nous cherchons les paramètres qui rendent nos observations les plus probables.

Construction de la vraisemblance¶

Supposons que nous avons un modèle paramétrique $p(y|\mathbf{x}; \boldsymbol{\theta})$ qui, pour chaque entrée $\mathbf{x}$ et choix de paramètres $\boldsymbol{\theta}$, définit une distribution sur les sorties possibles $y$. Par exemple, en régression, ce pourrait être une gaussienne centrée sur $f(\mathbf{x}; \boldsymbol{\theta})$.

Considérons un seul exemple $(\mathbf{x}_1, y_1)$. Pour des paramètres $\boldsymbol{\theta}$ fixés, nous pouvons évaluer $p(y_1 | \mathbf{x}_1; \boldsymbol{\theta})$: la probabilité (ou densité) que le modèle assigne à l’observation $y_1$. Si cette valeur est élevée, les paramètres $\boldsymbol{\theta}$ “expliquent bien” cette observation. Si elle est faible, $y_1$ est une valeur improbable sous ce modèle.